{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from catboost import Pool, CatBoostRegressor\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that we have two categorical features with 9 values each, so there are 81 possible feature combinations. The target y depends on the features x as y = mean(x) + eps(x), where mean(x) is some unknown fixed value and eps(x) is a normally distributed noise with mean 0 and variance var(x). We let mean(x) to be randomly generated and var(x) to have two values - 0.01 and 0.04. Also, there will be an \"unknown\" region without any training examples. To indicate this region on the figure, we set var(x) = 0 for these feature combinations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function for generating mean and variance\n",
    "def gen_parameters(noise=0.01, seed=0):\n",
    "    \n",
    "    np.random.seed(seed)\n",
    "    \n",
    "    mean = np.random.rand(9, 9)\n",
    "    \n",
    "    figure = np.array([\n",
    "        [1, 1, 1, 1, 1, 1, 1, 1, 1],\n",
    "        [1, 1, 4, 1, 1, 1, 4, 1, 1],\n",
    "        [1, 4, 0, 4, 1, 4, 0, 4, 1],\n",
    "        [4, 0, 0, 0, 4, 0, 0, 0, 4],\n",
    "        [1, 4, 0, 0, 0, 0, 0, 4, 1],\n",
    "        [1, 1, 4, 0, 0, 0, 4, 1, 1],\n",
    "        [1, 1, 1, 4, 0, 4, 1, 1, 1],\n",
    "        [1, 1, 1, 1, 4, 1, 1, 1, 1],\n",
    "        [1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
    "    ])\n",
    "\n",
    "    var = figure*noise\n",
    "    \n",
    "    return mean, var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True data uncertainty (variance). The white region will not contain any traing data.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAD8CAYAAACihcXDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAW7klEQVR4nO3df5BdZX3H8fdnNxAFNDjUOpggSUvUhtpBoUGrdaQpNqg1OoUxOK3USV2dmqp12orVItIfIx0KtQO1LoJFag2alnari9gB2mqLcQNCIUHqEkE2/uJXg4gIIZ/+cQ96vd299+7u/XHO2c9r5gznnh/3+9zZ5btPvud5nivbRETEcI0MuwEREZFkHBFRCknGERElkGQcEVECScYRESWQZBwRUQJJxhERc5C0UdLtkqYlnTnL+eWSrijO75C0uji+WtL3Jd1UbH/TKdayLhrzXGATsLI4tBeYsH3bfD5URESVSBoFLgJOBmaAKUkTtnc3XbYFeMD2MZI2A+cCryvO3WH7uG7jte0ZS3oXsA0Q8KViE/CJ2f5KRETUyHpg2vYe24/SyIWbWq7ZBFxW7G8HNkjSQoJ16hlvAY61/VjzQUnnA7uAD8x2k6QxYAzgQ+/eevzYa09ZSNsiYokZ+flXLCiRNXuLntr1tOIP8903U+Sqwrjt8WJ/JXB307kZ4MSWt/jhNbb3S9oHHFGcWyPpy8CDwHttf75dWzol4wPAM4G7Wo4fWZybVfFhxgEOTE1mvnVElFJzruqxbwLPsn2fpOOBf5J0rO0H57qhUzJ+B3CNpK/yo78QzwKOAbb2osUREb3Uw1EJe4Gjml6vKo7Nds2MpGXACuA+Nxb9+QGA7Rsk3QE8G9g5V7C2ydj2ZyU9m0btpPkB3pTtx7v+SBERA7JsYSXb2UwBayWtoZH3NgOvb7lmAjgDuB44FbjWtiU9Hbjf9uOSfgpYC+xp2+5OrbF9APjivD9GRMQQjPQoFxc14K3A1cAocKntXZLOAXbangAuAS6XNA3cTyNhA7wUOEfSYzRKum+xfX+7eOr3EpqpGUdEt3rxAO+doyu6zjnnP76vZ93oxerYM46IqJKR3pUpBirJOCJqparTipOMI6JWelUzHrQk44ioldGUKSIihi9lioiIEkiZIiKiBNIzDq57xZuGEvekyYuHEncpyc+2OjK0LSKiBJZVMxcnGUdEvaRMERFRAiNUs2ucZBwRtZLRFBERJZAyRURECaRnHBFRAj1cXH6gkowjolZSpoiIKIGUKSIiSiBD2yIiSiA944iIEhhNMo6IGL6qlikW/OBR0ht72ZCIiF4YUfdbmSxmFMj75zohaUzSTkk7x6+8ahEhIiLmZ2QeW5m0LVNI+u+5TgHPmOs+2+PAOMCBqUkvuHUREfNUsg5v1zrVjJ8B/ArwQMtxAf/VlxZFRCxCXReX/zRwmO2bWk9I+re+tCgiYhHKVn7oVttkbHtLm3Ov731zIiIWp5r94gxti4iaUU3LFBERlVLNVJxkHBE1U8uacURE1VS0SpFkHBH1UtXp0EnGEVEr1UzFScYRUTNlW3OiW0nGEVErqmjfuJbJ+LpXvGkocTfctXsoca85et1Q4p40efFQ4g7j55ufbXX0MhVL2gh8EBgFPmL7Ay3nlwMfA44H7gNeZ/vOpvPPAnYDZ9s+r12sqo4CiYiYVa+W0JQ0ClwEnAKsA06X1PrXcQvwgO1jgAuAc1vOnw90tXRlknFE1MoI6nrrYD0wbXuP7UeBbcCmlms2AZcV+9uBDSqmAEp6DfA1YFd37Y6IqBHNZ2tae73YxpreaiVwd9PrmeIYs11jez+wDzhC0mHAu2iz7nurWtaMI2Lpms+kj+a113vsbOAC2w91u1ZGknFE1EoPH+DtBY5qer2qODbbNTOSlgEraDzIOxE4VdKfA4cDByQ9YvvCuYIlGUdErfRwaNsUsFbSGhpJdzPQunTwBHAGcD1wKnCtbQO/+MP2SGcDD7VLxJBkHBE1M9qjXGx7v6StwNU0hrZdanuXpHOAnbYngEuAyyVNA/fTSNgLkmQcEbXSy3HGtieByZZjZzXtPwKc1uE9zu4mVpJxRNRKZuBFRJRAltCMiCiBqk6e6NhuSc+VtKEYxNx8fGP/mhURsTDzmfRRJm2TsaS3Af8M/A5wq6TmqYB/1s+GRUQsxIjU9VYmnXrGbwKOt/0a4GXAH0l6e3Fuzk/SPMVw/Mqu1siIiOiJqvaMO9WMR2w/BGD7TkkvA7ZLOpo2n6V5iuGBqUn3qK0RER11O/24bDr1jL8t6bgnXhSJ+VXATwDP62fDIiIWoldLaA5ap57xG4D9zQeKlYneIOnDfWtVRMQCqWxZtkttk7HtmTbn/rP3zYmIWJyRio5tyzjjiKiVqtaMk4wjolYqmouTjCOiXtIzjogogYrm4iTjiKiXss2s61aScUTUykgdh7ZFRFSNMrQtImL4qvoAT43vzuufa56+cuBrU2y4a/egQy5J1xy9bihx8/Ptv6H9bO/Zu+hMese6tV3nnJ/e/dXSZO70jCOiVqraM04yjohaqWguTjKOiHoZzWiKiIjhS5kiIqIEKpqLk4wjol6SjCMiSqCWi8tHRFRNHuBFRJRAyhQRESWQ0RQRESVQ0VzcORlLWg/Y9pSkdcBG4Cu2J/veuoiIeaplz1jS+4BTgGWS/hU4EbgOOFPS823/6QDaGBHRtYrm4o4941OB44DlwLeAVbYflHQesAOYNRlLGgPGAN5x2Ape9aRDe9fiiIg2RkarmY07JeP9th8HHpZ0h+0HAWx/X9KBuW6yPQ6Mw3CW0IyIpauqZYpOa+I/KumQYv/4Jw5KWgHMmYwjIoZmRN1vHUjaKOl2SdOSzpzl/HJJVxTnd0haXRxfL+mmYrtZ0ms7NrvD+ZfafhjAdnPyPQg4o+MniYgYNKn7re3baBS4iMZzs3XA6cUghmZbgAdsHwNcAJxbHL8VOMH2cTQGPXxYUttKRNtkbPsHcxy/1/YtbT9JRMQQSOp662A9MG17j+1HgW3AppZrNgGXFfvbgQ2SZPth2/uL408COpZrK/rVfRERcxgd6XqTNCZpZ9M21vROK4G7m17PFMeY7Zoi+e4DjgCQdKKkXcAtwFuakvOsMukjImplPgsFNQ826DXbO4BjJf0McJmkq2w/Mtf16RlHRL30qGYM7AWOanq9qjg26zVFTXgFcF/zBbZvAx4CfrZdsCTjiKgVjajrrYMpYK2kNZIOBjYDEy3XTPCjwQynAtfadnHPMgBJRwPPBe5sFyxlioiolx6NM7a9X9JW4GpgFLjU9i5J5wA7bU8AlwCXS5oG7qeRsAFeQmOm8mM0hgH/tu1728VLMo6IeunhesbFGjyTLcfOatp/BDhtlvsuBy6fT6wk44ioFY1Ws/qaZBwR9VLR6dB9T8YnTV7c7xD/zzVHt06SGYwNd+0eStxhWWqfdxiG9bs8jP9ve0XV7BinZxwRNZOecUTE8OXboSMiyiA944iI4ctoioiIMkiZIiKiBFKmiIgYvqp+7VKScUTUS8oUERHDlwd4ERFlkDJFRMTwVXXSx7z785I+1o+GRET0RO++6WOg2vaMJbWuai/gJEmHA9h+db8aFhGxIDXtGa8CHgTOB/6i2L7btD+r5m9cHb/yql61NSKiI0ldb2XSqWZ8AvB24D3A79u+SdL3bf97u5uav3H1wNSke9LSiIhu1HE0he0DwAWSPlX899ud7omIGKqS9Xi71VVitT0DnCbplTTKFhER5VTnZPwE258BPtOntkRELN5IDcsUERGVsxR6xhERpZdkHBFRAqOjw27BgiQZR0S9pGccEVECScYRESWQZBwRUQIZ2hYRUQJJxuVx0uTFQ4l7zdHrhhJ3w127hxJ3KRnWz3ZYv8uVljJFRMTwKT3jiIgSSM84IqIEkowjIkogyTgiogQqOh26mpXuiIi59PALSSVtlHS7pGlJZ85yfrmkK4rzOyStLo6fLOkGSbcU//2lTrGSjCOiXnqUjCWNAhcBpwDrgNMltY5x3AI8YPsY4ALg3OL4vcCv2n4ecAZweadmJxlHRL2MjHS/tbcemLa9x/ajwDZgU8s1m4DLiv3twAZJsv1l298oju8Cnixpedtmz+tDRkSU3Tx6xs3fZF9sY03vtBK4u+n1THGM2a6xvR/YBxzRcs2vATfa/kG7ZucBXkTUyzxGUzR/k31/mqJjaZQuXt7p2iTjiKiX3o2m2Asc1fR6VXFstmtmJC0DVgD3AUhaBVwJvMH2HZ2CzatMIeklkt4pqWOWj4gYit6NppgC1kpaI+lgYDMw0XLNBI0HdACnAtfatqTDaXx585m2/7ObZrdNxpK+1LT/JuBC4CnA+2Yb5hERMXQ9SsZFDXgrcDVwG/BJ27sknSPp1cVllwBHSJoG3gk8kRe3AscAZ0m6qdh+sl28TmWKg5r2x4CTbd8j6Tzgi8AHZrupKIKPAXzo3VsZe+0pHcJERPRIDxcKsj0JTLYcO6tp/xHgtFnu+xPgT+YTq1MyHpH0NBo9aNm+pwj0PUn757qpuSh+YGrS82lQRMSi1HQ69ArgBkCAJR1p+5uSDiuORUSUy0g1p0O3Tca2V89x6gDw2p63JiJisUaq2U9c0NA22w8DX+txWyIiFk/VnMuWccYRUS81rRlHRFRLvnYpIqIE0jOOiCiBOo6miIionJQpIiJKIGWKiIgSyNC2iIgSWEqTPmJ2J01ePJS41xzd+rVcg7Hhrt1DiTuMzzusn20sQB7gRUSUQMoUERElkDJFREQJZDRFREQJpEwREVECKVNERJRARlNERJRAyhQRESWQMkVERAmkZxwRUQIZ2hYRUQJZQjMiogQqOpqi7Z8QSSdKemqx/2RJ75f0L5LOlbRiME2MiJgHqfutRDr15y8FHi72PwisAM4tjn10rpskjUnaKWnn+JVX9aShERFdGRnpfiuRTmWKEdv7i/0TbL+g2P+CpJvmusn2ODAOcGBq0otvZkREl0rW4+1Wpz8Nt0p6Y7F/s6QTACQ9G3isry2LiFgIjXS/lUinnvFvAR+U9F7gXuB6SXcDdxfnIiLKpaIP8NomY9v7gN8sHuKtKa6fsf3tQTQuImLe6jwDz/aDwM19bktExOKVrPzQrYwzjoh6qegDvCTjiKiXivaMq9nqiIg5SOp66+K9Nkq6XdK0pDNnOb9c0hXF+R2SVhfHj5B0naSHJF3YTbuTjCOiXkaWdb+1IWkUuAg4BVgHnC5pXctlW4AHbB8DXEBjUhzAI8AfAb/XdbO7vTAiohJG1P3W3npg2vYe248C24BNLddsAi4r9rcDGyTJ9vdsf4FGUu6u2d1eGBFRCfOY9NG8dEOxjTW900oacyqeMFMcY7ZritnK+4AjFtLsPMCLiHqZx2iK5qUbhi3JOCLqpXejKfYCRzW9XlUcm+2aGUnLaCymdt9CgiUZ18BJkxcPJe41R7c+yxiMYX3eqIjejTOeAtZKWkMj6W4GXt9yzQRwBnA9cCpwre0FLY6WZBwR9TLam7UpbO+XtBW4GhgFLrW9S9I5wE7bE8AlwOWSpoH7aSRsACTdCTwVOFjSa4CX2949V7wk44iolx5O+rA9CUy2HDuraf8R4LQ57l09n1hJxhFRL5kOHRFRAhWdDp1kHBH1kp5xREQJjFYzrVWz1RERc+hmAaAySjKOiHpJzTgiogTSM46IKIH0jCMiSqCiPeO2f0IkvU3SUe2uiYgoldHR7rcS6dSf/2Ngh6TPS/ptSU8fRKMiIhZsHusZl0mn1uyhsWzcHwPHA7slfVbSGZKeMtdNzQs2j195VQ+bGxHRgdT9ViKdasa2fQD4HPA5SQfR+D6o04HzgFl7ys0LNh+YmlzQcnIREQtTriTbrU7J+Mc+le3HaKzfOSHpkL61KiJioUrW4+1Wp2T8urlO2H64x22JiFi8OiZj2/8zqIZERPREyR7MdSvjjCOiXqrZMU4yjoi6qWY2TjKOiHqpY804IqJykowjIkogD/AiIsogPeOIiOFLmSIiogSSjGOpOWny4mE3IWIWScYREUOXLySNiCiDjKaIiCiB9IwjIkogyTgiogySjCMihi8944iIEqhmLk4yjoiayWiKiIgSSJkiIqIMqpmMq9mfj4iYi9T91vGttFHS7ZKmJZ05y/nlkq4ozu+QtLrp3LuL47dL+pVOsdomY0kHS3qDpF8uXr9e0oWS3irpoI6fJCJi0HqUjCWNAhcBpwDrgNMlrWu5bAvwgO1jgAuAc4t71wGbgWOBjcBfF+83p05lio8W1xwi6QzgMOAfgQ3AeuCMDvdHRAxW7x7grQembe8BkLQN2ATsbrpmE3B2sb8duFCNxTE2Adts/wD4mqTp4v2unytYp2T8PNs/J2kZsBd4pu3HJf0dcPNcN0kaA8aKl2+2Pd4hzpzvs9B7F2MpxV1Kn3WpxV1Kn/XHHLKi66JxS64CGG9q+0rg7qZzM8CJLW/xw2ts75e0DziiOP7FlntXtmtLpz8hI5IOBp4CHAKsKI4vB+YsU9get31CsS3mhzLW+ZK+WEpxl9JnXWpxl9JnXZCWXLXYfLUonXrGlwBfAUaB9wCfkrQHeCGwrc9ti4gYpr3AUU2vVxXHZrtmpqggrADu6/LeH9O2Z2z7AuAlwIts/xXwa8DVwBbb7+/4USIiqmsKWCtpTVEh2AxMtFwzwY+enZ0KXGvbxfHNxWiLNcBa4EvtgnUcZ2z7G037/0ujSD0ow/onw1KKu5Q+61KLu5Q+a88VNeCtNDqgo8CltndJOgfYaXuCRvXg8uIB3f00EjbFdZ+k8bBvP/BW24+3i6dGEo+IiGHKpI+IiBJIMo6IKIHSJuNO0xD7FPNSSd+RdOsg4hUxj5J0naTdknZJevuA4j5J0pck3VzEHdgDWUmjkr4s6dODilnEvVPSLZJukrRzQDEPl7Rd0lck3SbpRQOI+ZziMz6xPSjpHf2OW8T+3eL36VZJn5D0pEHErYNS1oyLaYP/A5xMY7D0FHC67d1tb1x83JcCDwEfs/2z/YzVFPNI4EjbN0p6CnAD8JoBfFYBh9p+qJja/gXg7ba/2OHWXsR+J3AC8FTbr+p3vKa4dwIn2L53gDEvAz5v+yPFE/lDigfhg4o/SmNI1Ym27+pzrJU0fo/W2f5+8QBr0vbf9jNuXZS1Z/zDaYi2H6UxpnlTv4Pa/g8aT0QHxvY3bd9Y7H8XuI0OM3V6FNe2HypeHlRsff/LLGkV8ErgI/2ONWySVgAvpfHEHduPDjIRFzYAd/Q7ETdZBjy5GHN7CPCNDtdHoazJeLZpiH1PUMNWrPj0fGDHgOKNSroJ+A7wr7YHEfcvgT8ADgwgVisDn5N0QzENtt/WAPcAHy3KMh+RdOgA4jbbDHxiEIFs7wXOA74OfBPYZ/tzg4hdB2VNxkuOpMOAfwDeYfvBQcS0/bjt42jMDlovqa+lGUmvAr5j+4Z+xmnjJbZfQGMVrrcWZal+Wga8APiQ7ecD3wMG8vwDGqsuAq8GPjWgeE+j8S/YNcAzgUMl/fogYtdBWZPxvKcSVllRs/0H4OO2/3HQ8Yt/Ol9HY6m/fnox8OqidrsN+KVi0amBKHpu2P4OcCWNclg/zQAzTf/i2E4jOQ/KKcCNtr89oHi/DHzN9j22H6OxwuMvDCh25ZU1GXczDbEWigdplwC32T5/gHGfLunwYv/JNB6WfqWfMW2/2/Yq26tp/EyvtT2QnpOkQ4sHpBSlgpcDfR01Y/tbwN2SnlMc2sCPL7/Yb6czoBJF4evACyUdUvxeb6DxDCS6UMqvXZprGmK/40r6BPAy4CckzQDvs31Jn8O+GPgN4Jaifgvwh7Yn+xz3SOCy4mn7CPBJ2wMdajZgzwCubOQIlgF/b/uzA4j7O8DHi07FHuCNA4j5xB+ck4E3DyIegO0dkrYDN9KYAvxlajI1ehBKObQtImKpKWuZIiJiSUkyjogogSTjiIgSSDKOiCiBJOOIiBJIMo6IKIEk44iIEvg/4dRQO7diVgEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate parameters and plot variance\n",
    "mean, var = gen_parameters()\n",
    "print(\"True data uncertainty (variance). The white region will not contain any traing data.\")\n",
    "sns.heatmap(var, cmap=\"Reds\", vmax=0.05)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function for generating train and validation sets\n",
    "def generate_training_data(n_samples, mean, var, num_cat=9, seed=0):\n",
    "    np.random.seed(seed)\n",
    "    train = []\n",
    "    target = []\n",
    "    val = []\n",
    "    val_target = []\n",
    "    for i in range(num_cat):\n",
    "        for j in range(num_cat):\n",
    "            if var[i, j] == 0:\n",
    "                continue\n",
    "            for _ in range(n_samples):\n",
    "                train.append([i, j])\n",
    "                val.append([i, j])\n",
    "                target.append(np.random.normal(mean[i, j], math.sqrt(var[i, j])))\n",
    "                val_target.append(np.random.normal(mean[i, j], math.sqrt(var[i, j])))\n",
    "                \n",
    "    train = np.asarray(train)\n",
    "    target = np.asarray(target)\n",
    "\n",
    "    return train, target, val, val_target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we generate the dataset using the obtained distribution. In each cell, except for the white ones, we have several training examples. The test set consists of all possible feature combinations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate train and validation datasets\n",
    "train, target, val, val_target = generate_training_data(1000, mean, var)\n",
    "\n",
    "train_pool = Pool(train, target, cat_features = [0, 1])\n",
    "val_pool = Pool(val, val_target, cat_features = [0, 1])\n",
    "\n",
    "# generate test, consisting of all possible feature combinations\n",
    "num_cat = 9\n",
    "test = np.asarray([[i, j] for i in range(num_cat) for j in range(num_cat)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we optimize the RMSE loss, then we predict mean(x). If we also want to estimate var(x) (data uncertainty), then we have to use a probabilistic regression model that predicts mean and variance. For this purpose, we can use the loss function called RMSEWithUncertainty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best iteration = 978\n"
     ]
    }
   ],
   "source": [
    "# predict mean value and data uncertainty\n",
    "\n",
    "model = CatBoostRegressor(iterations=1000, learning_rate=0.2, loss_function='RMSEWithUncertainty',\n",
    "                          verbose=False, random_seed=0)\n",
    "\n",
    "model.fit(train_pool, eval_set=val_pool)\n",
    "print(\"best iteration =\", model.get_best_iteration())\n",
    "preds = model.predict(test)\n",
    "\n",
    "mean_preds = preds[:, 0] #the first prediction is the estimated mean value\n",
    "var_preds = preds[:, 1] #the second prediction is the estimated variance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see how well we can predict the mean values with this loss. Below we compare the true mean values with the estimated ones. For easier visual comparison, we use a mask to hide all feature combinations for which we do not have any training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define mask to hide elements inside the heart:\n",
    "mask = np.array([\n",
    "        [1, 1, 1, 1, 1, 1, 1, 1, 1],\n",
    "        [1, 1, 1, 1, 1, 1, 1, 1, 1],\n",
    "        [1, 1, 0, 1, 1, 1, 0, 1, 1],\n",
    "        [1, 0, 0, 0, 1, 0, 0, 0, 1],\n",
    "        [1, 1, 0, 0, 0, 0, 0, 1, 1],\n",
    "        [1, 1, 1, 0, 0, 0, 1, 1, 1],\n",
    "        [1, 1, 1, 1, 0, 1, 1, 1, 1],\n",
    "        [1, 1, 1, 1, 1, 1, 1, 1, 1],\n",
    "        [1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
    "    ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True mean values (with mask):\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD8CAYAAADUv3dIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAVqElEQVR4nO3df5RcZX3H8fdnJ4D8DFpEMIkQD0GN2CNhG6RwEMoPgyKxVluCVvBQ154aRK1arBYFWi0URU+l6gr4q0qEiBg1Clb5pYaQ5YeQhB+GEMgCElRCilBDkm//mBvOsGTnzuzOPPfO3c+Lc09m7r1zv88c9nz32e99nucqIjAzszT6im6AmdlE4qRrZpaQk66ZWUJOumZmCTnpmpkl5KRrZpaQk66ZWUKT8k6Q9HJgLjAl2/UgsCgi7uxmw8zMqqhpT1fSPwELAAE3ZZuASyWd0f3mmZlVi5rNSJN0D/DKiHh6xP7tgRURMWOUzw0AAwD/edB+B5360r061+IWrLppOGm8rV4yfXIhcbffc7di4r5yv0Lirlu0JHnMFx7XnzwmwMWf/VEhcX/91KZC4p6/ab3Ge42/124tT7P9YmwYd7x25dV0twAv3sb+vbNj2xQRgxHRHxH9qROumVmZ5dV03wf8VNKvgbXZvpcA+wHzu9kwM7OxKPvogKZJNyJ+LGl/YDbPvpG2LCI2d7txZmbtmqTkFYO25I5eiIgtwI0J2mJmNm595c65+UnXzKyX9HR5wcys1/T1ennBzKyXuKdrZpaQa7pmZgnVXF4wM0vH5QUzs4RcXjAzS2jC93S3P2RWt0M8x8xDZnHNJxcmj7v/tD9JHhNg0/onC4l73qeuLCTuGfffkj7olmIWgHn3hy8oJO6V+xxQSNxO8JCxAhSRcM2sHCaVO+dWM+ma2cQ14csLZmYp9VHurq6TrplVikcvmJkl5PKCmVlC7umamSXU84uYm5n1EpcXzMwScnnBzCwhDxkzM0vIPV0zs4RqTrpmZumUvbww5ht9kt7ZyYaYmXVCn1rfCmnfOD571mgHJA1IGpI09OUly8cRwsysPX1tbEVoWl6QdPtoh4AXjfa5iBgEBgE2X3B6jLl1ZmZtKndxIb+m+yLgdcBjI/YL+GVXWmRmNg69voj5D4BdIuK2kQckXduVFpmZjUNPz0iLiFObHDup880xMxufcvdzPWTMzCpGPV5eMDPrKeVOuU66ZlYxPV3TNTPrNSWvLjjpmlm1lH0asJOumVVKuVNu+csfZmZt6eTaC5LmSLpb0ipJZ2zj+EskXSPpVkm3S3p9bvvG9rXMzMpJbfzX9DpSDbgQOA6YCcyTNHPEaR8DLouIA4ETgf/Ka1/Xyws3fOo73Q7xHA9tfDp5TIAdvnhFIXGv2+cVhcT90GlHFxL32gMOTR7zyDUrk8cEWHfEIYXEnbv0+4XE7YQOlhdmA6siYjWApAXAXKDxhyGA3bLXk4GH8i7qmq6ZVUo7SzZKGgAGGnYNZgt2AUwB1jYcGwYOHnGJTwBXSzoN2BnI7Yk46ZpZpbQzeqFxRcQxmgd8NSI+LekQ4BuSDoiILaO3z8ysQtTGluNBYFrD+6nZvkanApcBRMQS4HnAHs0u6qRrZpUitb7lWAbMkDRd0vbUb5QtGnHOA8BR9bh6BfWk+2izi7q8YGaV0qkbaRGxSdJ84CqgBlwSESsknQ0MRcQi4B+BL0t6P/WbaqdERNMHNzjpmlml5A0Fa0dELAYWj9h3ZsPrlUBbw2mcdM2sUvwIdjOzhEqec510zaxaOlle6AYnXTOrFC/taGaWUNnHwea2T9LLJR0laZcR++d0r1lmZmPTwckRXdE06Up6L/A94DRguaS5DYc/2c2GmZmNRZ/U8lZI+3KOvws4KCLeBBwB/Iuk07Njo7ZY0oCkIUlD33/qD51pqZlZC8re082r6fZFxBMAEbFG0hHAQkn70KTNjYtIXLvn1KazM8zMOqnsj2DP6+k+IunVW99kCfh46gs6vKqbDTMzG4tOPjmiG/J6uu8ANjXuiIhNwDskfalrrTIzGyMVlU1b1DTpRsRwk2O/6HxzzMzGp6/kY8Y8TtfMKqXsNV0nXTOrlJLnXCddM6sW93TNzBIqec510jWzailqplmrnHTNrFL6ennImJlZr5GHjJmZpTPhb6T9dH36BW/OWb8mecwiHXv/nYXEvX76KwuJe+SalYXELcKe1y4pJO761x1WSNzdb7hj3Ncoec51T9fMqmXC93TNzFIqec510jWzaql59IKZWTouL5iZJVTynOuka2bV4qRrZpZQTy9ibmbWa3wjzcwsIZcXzMwS8ugFM7OESp5z85OupNlARMQySTOBOcBdEbG4660zM2tTT/d0JX0cOA6YJOknwMHANcAZkg6MiH9L0EYzs5aVPOfm9nTfArwa2AH4DTA1IjZIOh9YCmwz6UoaAAYAXl/bkVl9O3SuxWZmTfTVyp1185b73RQRmyPiSeDeiNgAEBFPAVtG+1BEDEZEf0T0O+GaWUqSWt5auNYcSXdLWiXpjFHO+WtJKyWtkPStvGvm9XQ3StopS7oHNQSZTJOka2ZWmA6N05VUAy4EjgGGgWWSFkXEyoZzZgAfAQ6NiMck7ZnbvJzjh2cJl4hoTLLbASe3+R3MzLpPan1rbjawKiJWR8RGYAEwd8Q57wIujIjHACJiXd5FmybdiPjjKPt/GxHjX+LdzKzD2ikvSBqQNNSwDTRcagqwtuH9cLav0f7A/pJ+IelGSXPy2udxumZWLbXWn0wZEYPA4DiiTQJmAEcAU4HrJb0qItY3+4CZWWV0cMGbB4FpDe+nZvsaDQNLI+Jp4D5J91BPwstGu2jJH1ZsZtamztV0lwEzJE2XtD1wIrBoxDlXUu/lImkP6uWG1c0u6p6umVVKp3q6EbFJ0nzgKqAGXBIRKySdDQxFxKLs2LGSVgKbgQ9FxO+aXddJ18yqpYNT0rLlDhaP2Hdmw+sAPpBtLXHSNbNq8Xq6ZmbpqI3RC0Vw0jWzain5ijddT7offM20/JM6bP3rDkseE2D3q35eSNyiHH7fiqKbUH1Pbigk7IbHNxYSd/cOXEPl7ui6p2tmFTPRe7pmZin5acBmZim5p2tmlo5HL5iZpeTygplZQi4vmJml09NPAzYz6zkuL5iZpeMbaWZmKbm8YGaWTtknR7TdD5f09W40xMysIzr35IiuaNrTlTTy0RQCjpS0O0BEnNCthpmZjUnJe7p55YWpwErgIiCoJ91+4NPNPpQ9xngA4LP7TeGUvV8w/paambWg7EPG8soL/cDNwEeBxyPiWuCpiLguIq4b7UMRMRgR/RHR74RrZknV+lrfCtC0pxsRW4ALJF2e/ftI3mfMzApV8p5uSwk0IoaBt0p6A1DMqspmZq2oQtLdKiJ+CPywS20xMxu/Pk+OMDNLp0o9XTOz0nPSNTNLqFYrugVNOemaWbW4p2tmlpCTrplZQk66ZmYJeciYmVlCEz3prn/sj90O8RzTBs9LHhMg1q0pJK723LeQuBPJpk/NLyTu/H/9XiFxvzC8rJC4HeHygplZOproPV0zs6Tc0zUzS8hJ18wsoZIn3XIXP8zM2lWrtb7lkDRH0t2SVkk6o8l5fyUpJPXnXdNJ18yqpUMPppRUAy4EjgNmAvMkzdzGebsCpwNLW2mek66ZVUvnngY8G1gVEasjYiOwAJi7jfPOAc4F/q+V5jnpmlm19PW1vEkakDTUsA00XGkKsLbh/XC27xmSZgHTsgc8tMQ30sysWtq4kRYRg8Dg2MKoD/gMcEo7n3PSNbNq6dzohQeBaQ3vp2b7ttoVOAC4Nnvs+17AIkknRMTQaBd10jWzauncIubLgBmSplNPticCJ209GBGPA3tsfS/pWuCDzRIutJl0JR1Gvbi8PCKubuezZmZJdKinGxGbJM0HrgJqwCURsULS2cBQRCway3WbJl1JN0XE7Oz1u4D3AN8FPi5pVkT8+1iCmpl1TQcnR0TEYmDxiH1njnLuEa1cM2/0wnYNrweAYyLiLOBY4G2jfajxjuC3fr++lXaYmXVGG6MXipBXXuiT9HzqyVkR8ShARPxB0qbRPtR4R/D+V70sOtVYM7NcJZ8GnJd0JwM3AwJC0t4R8bCkXbJ9Zmbl0tfDTwOOiH1HObQF+MuOt8bMbLz6yt0fHNOQsYh4Erivw20xMxs/lXuircfpmlm19HhN18yst/hxPWZmCbmna2aWUC+PXjAz6zkuL5iZJeTygplZQh4yZmaWUMknRyiiu0sjbFnyvfRrL0zeI/+cLnj4bacWEnevc0d9SGlX1Q4rZlLi5iVjWlFvXLT3S5PHBNBuLygk7kNvLOb/7ZRb7xp3xtx86Xkt55zavA8nz9Du6ZpZtbi8YGaWUMnLC066ZlYtHr1gZpaQywtmZgm5vGBmlpCnAZuZJeTygplZQi4vmJkl5J6umVlCHjJmZpaQl3Y0M0uo5KMXmv5KkHSwpN2y1ztKOkvS9yWdK2lymiaambVBan0rQF4//BLgyez154DJwLnZvq+M9iFJA5KGJA0NXnlVRxpqZtaSvr7WtwLklRf6ImJT9ro/ImZlr38u6bbRPhQRg8AgFLS0o5lNXCW/kZaX6pdLemf2+leS+gEk7Q883dWWmZmNhfpa3wqQ19P9O+Bzkj4G/BZYImktsDY7ZmZWLiW/kdY06UbE48Ap2c206dn5wxHxSIrGmZm1rQoz0iJiA/CrLrfFzGz8PCPNzCyhkt9Ic9I1s2opeU+33K0zM2uTpJa3Fq41R9LdklZJes5jtyV9QNJKSbdL+qmkffKu6aRrZtXSN6n1rQlJNeBC4DhgJjBP0swRp91KfQ7DnwILgfNymzemL2VmVlZ9an1rbjawKiJWR8RGYAEwt/GEiLgmIrbO2r0RmJrbvDF8JTOz8mpjckTjkgXZNtBwpSnU5yRsNZztG82pwI/ymucbaWZWLW2MXmhcsmB8IfV2oB94bd65TrpmVi2dG73wIDCt4f3UbN+zw0lHAx8FXhsRf8y7aNeTrvacln9Sh8VjxUyY2+vU4wuJu3jeRwqJe/wdhxQS98q3fih5zDcvvz55TIDTp/5ZIXE/8/6jC4nbEZ0bp7sMmCFpOvVkeyJw0rND6UDgS8CciFjXykXd0zWzaql1Zu2FiNgkaT5wFVADLomIFZLOBoYiYhHwH8AuwOXZELQHIuKEZtd10jWzaung5IiIWAwsHrHvzIbXbf9J4KRrZtXiacBmZgmVfBqwk66ZVYt7umZmCdXKndbK3Tozsza1spBNkZx0zaxaXNM1M0vIPV0zs4Tc0zUzS6jkPd2mvxIkvVdS+sUTzMzGqlZrfStAXj/8HGCppBsk/YOkF6ZolJnZmLWxnm4R8qKupr6c2TnAQcBKST+WdLKkXUf7UOPCwIMLruhgc83MckitbwXIq+lGRGwBrgaulrQd9ecFzQPOB7bZ821cGDjuvSU611wzszzlrunmJd1ntT4ingYWAYsk7dS1VpmZjVXJb6TlJd2/Ge1Aw8PYzMzKo5eTbkTck6ohZmYd4XG6ZmYJlbuj66RrZlVT7qzrpGtm1dLLNV0zs57jpGtmlpBvpJmZpeSerplZOi4vmJklNNGT7ukzjux2iOc457XTk8cE+PVdvy8k7htX3VpI3M3XLSwk7ptv+5/kMbf8cnHymACfW/PLQuJuuWFRIXE7Y4InXTOzlPxgSjOzlDx6wcwsIfd0zcwSctI1M0vJSdfMLB33dM3MEip3znXSNbOK8egFM7OEXF4wM0up3Em33P1wM7N2Sa1vuZfSHEl3S1ol6YxtHN9B0rez40sl7Zt3zaZJV9L2kt4h6ejs/UmSPi/pPZK2y22xmVlqHUq6kmrAhcBxwExgnqSZI047FXgsIvYDLgDOzWteXnnhK9k5O0k6GdgFuAI4CpgNnJwXwMwsqc7dSJsNrIqI1QCSFgBzgZUN58wFPpG9Xgh8XpIiIka9akSMugG3Z/9OAh4Batl7bT02yucGgKFsG2gWIyf+mD87nm0ixZ1I33WixZ1I33U8bW3IVc/KV8BbgIsa3v8t8PkRn18OTG14fy+wR7OYeb8S+iRtD+wK7ARMzvbvAIxaXoiIwYjoz7bBnBjNDIzjs+MxkeJOpO860eJOpO86JiNy1XjzVUvyygsXA3cBNeCjwOWSVgOvARZ0uW1mZkV6EJjW8H5qtm9b5wxLmkS9Y/q7ZhdtmnQj4gJJ385ePyTp68DRwJcj4qb22m9m1lOWATMkTaeeXE8EThpxziLq97aWUC9H/CyyOsNocsfpRsRDDa/XUy8Wp9L1rr7jTqjvOtHiTqTv2nERsUnSfOAq6n/tXxIRKySdDQxFxCLq1YBvSFoF/J56Ym5KOUnZzMw6yJMjzMwSctI1M0uotEk3b/pdl2JeImmdpOUp4mUxp0m6RtJKSSsknZ4o7vMk3STpV1ncs1LEzWLXJN0q6QepYmZx10i6Q9JtkoYSxdxd0kJJd0m6U9IhCWK+LPuOW7cNkt7X7bhZ7PdnP0/LJV0q6Xkp4vaSUtZ0s+l39wDHAMPU7yLOi4iVTT84/riHA08AX4+IA7oZqyHm3sDeEXGLpF2Bm4E3JfiuAnaOiCeyKd0/B06PiBu7GTeL/QGgH9gtIo7vdryGuGuA/oj4bcKYXwNuiIiLsjHvO2U3pFPFr1G/835wRNzf5VhTqP8czYyIpyRdBiyOiK92M26vKWtP95npdxGxkfqY4LndDhoR11O/A5lMRDwcEbdkr/8XuBOYkiBuRMQT2dvtsq3rv4ElTQXeAFzU7VhFkzQZOJz6HW4iYmPKhJs5Cri32wm3wSRgx2zM6k7AQznnTzhlTbpTgLUN74dJkIiKlq1QdCCwNFG8mqTbgHXATyIiRdzPAh8GtiSINVIAV0u6WVKKWVPTgUeBr2TllIsk7ZwgbqMTgUtTBIqIB4HzgQeAh4HHI+LqFLF7SVmT7oQjaRfgO8D7ImJDipgRsTkiXk19ps1sSV0tqUg6HlgXETd3M04Th0XELOqrRr0nKyd10yRgFvCFiDgQ+AOQ5P4E1FcJBE4ALk8U7/nU/yKdDrwY2FnS21PE7iVlTbqtTL+rjKym+h3gmxFxRer42Z+81wBzuhzqUOCErLa6APgLSf/d5ZjPyHpiRMQ64LvUy1jdNAwMN/wFsZB6Ek7lOOCWiHgkUbyjgfsi4tGIeJr6ioR/nih2zyhr0n1m+l322/pE6tPtKie7oXUxcGdEfCZh3BdK2j17vSP1m5Z3dTNmRHwkIqZGxL7U/5/+LCKS9IQk7ZzdqCT7E/9Y6itEdU1E/AZYK+ll2a6jePaygN02j0SlhcwDwGsk7ZT9XB9F/R6FNSjl43pGm37X7biSLgWOAPaQNAx8PCIu7nLYQ6kvGXdHVl8F+OeIWNzluHsDX8vubvcBl0VE0iFcib0I+G49FzAJ+FZE/DhB3NOAb2adh9XAOxPE3PqL5Rjg3SniAUTEUkkLgVuATcCtVGRKcCeVcsiYmVlVlbW8YGZWSU66ZmYJOemamSXkpGtmlpCTrplZQk66ZmYJOemamSX0/0A/HvQhIM58AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted mean values (with mask):\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD8CAYAAADUv3dIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAVnElEQVR4nO3de5RdZXnH8e9vJiByxyKCSYAsCWC8lMusIMLCKKChIrGt1YSq4EKnXRJEbcUoFgUUFxRFV6GWAfFWJQJCDRoBlwYR5JLhKgm3EAOZQAkoEBEkJHn6x9mxhyFz9jkz57x7nz2/D2svztm35x3IeubNs9/33YoIzMwsjZ6iG2BmNp446ZqZJeSka2aWkJOumVlCTrpmZgk56ZqZJeSka2aW0IS8EyTtDcwCJma7VgELIuKeTjbMzKyKGvZ0JX0amA8IuCXbBFwsaV7nm2dmVi1qNCNN0v3A6yLihWH7NweWRMTUEa7rB/oB/qNv6v7H7bFL+1rchHuvW5E03kav2esVhcSd8IqtCon7sr/eq5C4j/7ohuQxd37n/sljApz/9asLifvQn9cVEvesdU9prPf4Z23b9DTb/4o1Y47Xqrya7gbg1ZvYv0t2bJMiYiAi+iKiL3XCNTMrs7ya7seBX0h6AFiZ7dsV2AOY28mGmZmNRtlHBzRMuhFxlaQ9gem8+EHa4ohY3+nGmZm1aoKSVwxakjt6ISI2ADclaIuZ2Zj1lDvn5iddM7Nu0tXlBTOzbtPT7eUFM7Nu4p6umVlCrumamSXU6/KCmVk6Li+YmSXk8oKZWULjvqe7+fR9Oh3iJd44fR9+ccalyeO+duIOyWMCrH/m+ULinnXGFYXEnffw7emDvlDMf+PjT/pqIXGves2+hcRtBw8ZK0ARCdfMymFCuXNuNZOumY1f4768YGaWUg/l7uo66ZpZpXj0gplZQi4vmJkl5J6umVlCXb+IuZlZN3F5wcwsIZcXzMwS8pAxM7OE3NM1M0uo10nXzCydspcXRv2gT9KH2tkQM7N26FHzWyHtG8O1p450QFK/pEFJgxfcePcYQpiZtaanha0IDcsLku4a6RDwqpGui4gBYABg/TknxqhbZ2bWonIXF/Jruq8C3gE8OWy/gN90pEVmZmPQ7YuY/wTYOiLuGH5A0rUdaZGZ2Rh09Yy0iDiuwbGj298cM7OxKXc/10PGzKxi1OXlBTOzrlLulOuka2YV09U1XTOzblPy6oKTrplVS2WnAZuZlZFa2HLvJc2UdJ+kZZLmbeL4rpIWSbpd0l2S/ibvnk66ZlYp7Vp7QVIvcB5wBDANmCNp2rDTPgdcEhH7ArOB/8xt32h+KDOzslIL/+SYDiyLiOURsRaYD8wadk4A22aftwMeybtpx2u61335sk6HeIknXliXPCbAFudfUUjcKyfvXUjcT330rYXEXfS6NyeP+bYVS5PHBHjsLQcWEvcdN11ZSNx2aKWiK6kf6K/bNZCtHQMwEVhZd2wIOGDYLb4AXCPpBGAr4LC8mH6QZmaV0sqSjfWLc43SHODbEfEVSQcC35P0+ojYMNIFTrpmViltHL2wCphc931Stq/eccBMgIi4UdIWwI7A6pHbZ2ZWIW0cvbAYmCppiqTNqT0oWzDsnIeBQwEkvRbYAni80U3d0zWzSmnX5IiIWCdpLnA10AtcFBFLJJ0GDEbEAuBfgAskfYLaQ7VjI6LhGuJOumZWKe2cGhERC4GFw/adUvd5KXBQK/d00jWzSmliKFihnHTNrFL8CnYzs4RKnnOddM2sWlxeMDNLyEs7mpklVPbJB7ntk7S3pEMlbT1s/8zONcvMbHTaubRjJzRMupI+BvwYOAG4W1L9CjtndLJhZmaj0SM1vRXSvpzjHwH2j4h3AzOAf5N0YnZsxBZL6pc0KGnwJ8/9qT0tNTNrQtl7unk13Z6IeAYgIlZImgFcJmk3GrS5fuWeRTtNbDglzsysncr+Cva8nu5jkvbZ+CVLwEdSW0XnDZ1smJnZaLTrzRGdktfT/SDwohXBI2Id8EFJ53esVWZmo6SismmTGibdiBhqcOyG9jfHzGxseko+ZszjdM2sUspe03XSNbNKKXnOddI1s2pxT9fMLKGS51wnXTOrlqJmmjXLSdfMKqWnm4eMmZl1G3nImJlZOmV/kKactwWP2Wc32z752gtnPP1Q6pDj0rW7v66QuDNWLCkk7njy5OEHFxJ3hxt+O+aM+eC0qU3nnNcsfSB5hnZP18wqpew9XSddM6uUkudcJ10zq5Zej14wM0vH5QUzs4RKnnOddM2sWpx0zcwS6upFzM3Muo0fpJmZJeTygplZQh69YGaWUMlzbn7SlTQdiIhYLGkaMBO4NyIWdrx1ZmYt6uqerqTPA0cAEyT9HDgAWATMk7RvRHwpQRvNzJpW8pxL3sqT7wEOAg4BjgfeHRGnA+8A3jfSRZL6JQ1KGrx9w9q2NdbMLE9Pr5re8kiaKek+ScskzRvhnPdKWippiaQf5N0zr7ywLiLWA89KejAi1gBExHOSNox0UUQMAANQzNKOZjZ+tau8IKkXOA84HBgCFktaEBFL686ZCnwGOCginpS0U95983q6ayVtmX3evy7QdsCISdfMrDA9an5rbDqwLCKWR8RaYD4wa9g5HwHOi4gnASJidW7zco4fEhHPZjerT7KbAcfk3dzMLDmp6a2+FJpt/XV3mgisrPs+lO2rtyewp6QbJN0kaWZe8xqWFyLi+RH2PwE8kXdzM7PUWikv1JdCR2kCMBWYAUwCrpP0hoh4qtEFZmbV0du2N1OuAibXfZ+U7as3BNwcES8Av5N0P7UkvHikm5b8vZlmZq1Rj5reciwGpkqaImlzYDawYNg5/0Otl4ukHamVG5Y3uql7umZWLW0avRAR6yTNBa4GeoGLImKJpNOAwYhYkB17u6SlwHrgUxHx+0b3ddI1s0pp59KO2czbhcP2nVL3OYBPZltTnHTNrFpKPiXNSdfMqsXr6ZqZpaP2jV7oCCddM6uW8V5e+PSBu3Y6xEusnnFg8pgAO117YyFxizJjxZKim1B9z/2xkLBr1hSzUNUObbiHyt3RdU/XzCpmvPd0zcxS8tuAzcxSck/XzCwdj14wM0vJ5QUzs4RcXjAzS6er3wZsZtZ1XF4wM0vHD9LMzFJyecHMLJ2yT45ouR8u6budaIiZWVu08DbgIjTs6Uoa/j4gAW+VtD1ARBzVqYaZmY1KyXu6eeWFScBS4EIgqCXdPuArjS7K3h3fD/C1PSZy7C6vGHtLzcyaUPYhY3nlhT7gVuBk4OmIuBZ4LiJ+FRG/GumiiBiIiL6I6HPCNbOkenua3wrQsKcbERuAcyRdmv37sbxrzMwKVfKeblMJNCKGgH+Q9E5gTWebZGY2BlVIuhtFxE+Bn3aoLWZmY9fjyRFmZulUqadrZlZ6TrpmZgn19hbdgoacdM2sWtzTNTNLyEnXzCwhJ10zs4Q8ZMzMLKHxnnSfevL5Tod4ickXnJ08JkCsXlFIXO20eyFxx5N1X55bSNy5X/xxIXG/8fBNhcRtC5cXzMzS0Xjv6ZqZJeWerplZQiVPuuXuh5uZtaqNr+uRNFPSfZKWSZrX4Ly/lxSS+vLu6Z6umVVLm6YBS+oFzgMOB4aAxZIWRMTSYedtA5wI3NzMfd3TNbNqaV9PdzqwLCKWR8RaYD4waxPnnQ6cCfy5meY56ZpZtbSQdCX1Sxqs2/rr7jQRWFn3fSjbVxdK+wGTs7XGm+LygplVSwtDxiJiABgYTRhJPcBXgWNbuc5J18yqpX2jF1YBk+u+T8r2bbQN8Hrg2uwNxDsDCyQdFRGDI93USdfMqqV9SXcxMFXSFGrJdjZw9MaDEfE0sOP/h9W1wL82SrjgpGtmVdOm0QsRsU7SXOBqoBe4KCKWSDoNGIyIBaO5b0tJV9LB1J7o3R0R14wmoJlZR7VxckRELAQWDtt3ygjnzmjmng0rzpJuqfv8EeBcanWMzzcaKGxmVpg2To7ohLzHfJvVfe4HDo+IU4G3A/840kX1wzB+8Ien2tBMM7Mm9fQ0vxUgr7zQI2kHaslZEfE4QET8SdK6kS6qH4bx0Bv2inY11swsV8nXXshLutsBtwICQtIuEfGopK2zfWZm5dLTxW8DjojdRzi0AfjbtrfGzGysesrdHxzVkLGIeBb4XZvbYmY2dir36gYep2tm1dLlNV0zs+7i1/WYmSXknq6ZWULdPHrBzKzruLxgZpaQywtmZgl5yJiZWUIlnxyhiM4ujbDhN1ekX3th279KHhLg0Q/055/UATuffXIhcXsPPKqQuOt/fXnymNpt7+QxAbTVdoXEXfWuvysk7qQ77h1zxlx/8VlN55zeOSclz9Du6ZpZtbi8YGaWUMnLC066ZlYtHr1gZpaQywtmZgm5vGBmlpCnAZuZJeTygplZQi4vmJkl5J6umVlCHjJmZpaQl3Y0M0uo5KMXGv5KkHSApG2zzy+XdKqkKyWdKamYlTjMzBqRmt8KkNcPvwh4Nvv8dWA74Mxs37dGukhSv6RBSYMDP76mLQ01M2tKT0/zWwHyygs9EbEu+9wXEftln6+XdMdIF0XEADAABS3taGbjV8kfpOWl+rslfSj7fKekPgBJewIvdLRlZmajoZ7mtwLk9XQ/DHxd0ueAJ4AbJa0EVmbHzMzKpeQP0hom3Yh4Gjg2e5g2JTt/KCIeS9E4M7OWVWFGWkSsAe7scFvMzMbOM9LMzBLq8gdpZmbdpY0P0iTNlHSfpGWS5m3i+CclLZV0l6RfSNot755OumZWKZKa3nLu0wucBxwBTAPmSJo27LTbqQ2nfSNwGXBWXvucdM2sWnomNL81Nh1YFhHLI2ItMB+YVX9CRCyKiI0TyG4CJuXd1DVdM6uW9o1emEhteOxGQ8ABDc4/DvhZ3k2ddM2sWloYvSCpH+iv2zWQzahtLaT0fqAPeEveuU66ZlYtLYxeqF+yYBNWAZPrvk/K9g0Lp8OAk4G3RMTzeTGddM2sWto3TncxMFXSFGrJdjZw9ItCSfsC5wMzI2J1MzfteNLVTrt2OsRLxNNPJI8JsPOH31VI3IXvPamQuEfe2VdI3Eve9+nkMWffc0PymACfmFTMf+OzTzy0kLht0aZxuhGxTtJc4GqgF7goIpZIOg0YjIgFwL8DWwOXZqMhHo6Ioxrd1z1dM6uW3vatvRARC4GFw/adUvf5sFbv6aRrZtXiacBmZgmVfBqwk66ZVYt7umZmCbmna2aWUG+501q5W2dm1qK8hWyK5qRrZtXimq6ZWULu6ZqZJeSerplZQiXv6Tb8lSDpY5ImNzrHzKxUenub3wqQ1w8/HbhZ0q8lfVTSK1M0ysxs1Nr4jrROyIu6nNoakqcD+wNLJV0l6RhJ24x0kaR+SYOSBgfmX97G5pqZ5ZCa3wqQV9ONiNgAXANcI2kzai9pmwOcDWyy51u/MHAsuzXa11wzszzlrunmJd0XtT4iXgAWAAskbdmxVpmZjVbJH6TlJd33jXSg7g2YZmbl0c1JNyLuT9UQM7O28DhdM7OEyt3RddI1s6opd9Z10jWzaunmmq6ZWddx0jUzS8gP0szMUnJP18wsHZcXzMwSGu9J94Q939bpEC/xpYN3Sx4T4IEHniok7rvuHywk7vpbflZI3Nl3LUoec8P1VyaPCXDO8usLibvh+p8UErc9xnnSNTNLyS+mNDNLyaMXzMwSck/XzCwhJ10zs5ScdM3M0nFP18wsoXLnXCddM6sYj14wM0vI5QUzs5TKnXTL3Q83M2uV1PyWeyvNlHSfpGWS5m3i+Msk/TA7frOk3fPu2TDpStpc0gclHZZ9P1rSuZKOl7RZbovNzFJrU9KV1AucBxwBTAPmSJo27LTjgCcjYg/gHODMvObllRe+lZ2zpaRjgK2By4FDgenAMXkBzMySat+DtOnAsohYDiBpPjALWFp3zizgC9nny4BzJSkiYsS7RsSIG3BX9u8JwGNAb/ZdG4+NcF0/MJht/Y1i5MQf9bVj2cZT3PH0s463uOPpZx1LW+ty1YvyFfAe4MK67x8Azh12/d3ApLrvDwI7NoqZ9yuhR9LmwDbAlsB22f6XASOWFyJiICL6sm0gJ0Yj/WO4dizGU9zx9LOOt7jj6WcdlWG5aqz5qil55YVvAvcCvcDJwKWSlgNvAuZ3uG1mZkVaBUyu+z4p27epc4YkTaDWMf19o5s2TLoRcY6kH2afH5H0XeAw4IKIuKW19puZdZXFwFRJU6gl19nA0cPOWUDt2daN1MoRv4yszjCS3HG6EfFI3eenqBWLU+l4V99xx9XPOt7ijqefte0iYp2kucDV1P62f1FELJF0GjAYEQuoVQO+J2kZ8Adqibkh5SRlMzNrI0+OMDNLyEnXzCyh0ibdvOl3HYp5kaTVku5OES+LOVnSIklLJS2RdGKiuFtIukXSnVncU1PEzWL3SrpdUtJXzkpaIem3ku6QlOQVypK2l3SZpHsl3SPpwAQx98p+xo3bGkkf73TcLPYnsj9Pd0u6WNIWKeJ2k1LWdLPpd/cDhwND1J4izomIpQ0vHHvcQ4BngO9GxOs7Gasu5i7ALhFxm6RtgFuBdyf4WQVsFRHPZFO6rwdOjIibOhk3i/1JoA/YNiKO7HS8urgrgL6IeCJhzO8Av46IC7Mx71tmD6RTxe+l9uT9gIh4qMOxJlL7czQtIp6TdAmwMCK+3cm43aasPd2/TL+LiLXUxgTP6nTQiLiO2hPIZCLi0Yi4Lfv8R+AeYGKCuBERz2RfN8u2jv8GljQJeCdwYadjFU3SdsAh1J5wExFrUybczKHAg51OuHUmAC/PxqxuCTySc/64U9akOxFYWfd9iASJqGjZCkX7Ajcnitcr6Q5gNfDziEgR92vAScCGBLGGC+AaSbdKSjFragrwOPCtrJxyoaStEsStNxu4OEWgiFgFnA08DDwKPB0R16SI3U3KmnTHHUlbAz8CPh4Ra1LEjIj1EbEPtZk20yV1tKQi6UhgdUTc2sk4DRwcEftRWzXq+Kyc1EkTgP2Ab0TEvsCfgCTPJ6C2SiBwFHBpong7UPsb6RTg1cBWkt6fInY3KWvSbWb6XWVkNdUfAd+PiMtTx8/+yrsImNnhUAcBR2W11fnA2yT9d4dj/kXWEyMiVgNXUCtjddIQMFT3N4jLqCXhVI4AbouIxxLFOwz4XUQ8HhEvUFuR8M2JYneNsibdv0y/y35bz6Y23a5ysgda3wTuiYivJoz7SknbZ59fTu2h5b2djBkRn4mISRGxO7X/p7+MiCQ9IUlbZQ8qyf6K/3ZqK0R1TET8L7BS0l7ZrkN58bKAnTaHRKWFzMPAmyRtmf25PpTaMwqrU8rX9Yw0/a7TcSVdDMwAdpQ0BHw+Ir7Z4bAHUVsy7rdZfRXgsxGxsMNxdwG+kz3d7gEuiYikQ7gSexVwRS0XMAH4QURclSDuCcD3s87DcuBDCWJu/MVyOPBPKeIBRMTNki4DbgPWAbdTkSnB7VTKIWNmZlVV1vKCmVklOemamSXkpGtmlpCTrplZQk66ZmYJOemamSXkpGtmltD/AYQeGK7YMKt0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"True mean values (with mask):\")\n",
    "sns.heatmap(mean*mask, cmap=\"Reds\")\n",
    "plt.show()\n",
    "\n",
    "print(\"Predicted mean values (with mask):\")\n",
    "sns.heatmap(mean_preds.reshape([num_cat, num_cat])*mask, cmap=\"Reds\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that CatBoost captures the mean values well. Now let us plot the predicted data uncertainty (variance). We can see that CatBoost is able to predict the variance outside the heart and on its border. Inside the heart we have no training data, so anything can be predicted there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted data uncertainty:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAD8CAYAAACihcXDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYH0lEQVR4nO3df7RdZX3n8ffn3kAggqFDHQcTbDIDakNtUWJoR0TbFA2DY3SAMbBmpE7sxVVppa6ZKR1bRPpjZBYj4ywYNWOwSKtBg3SuNPywRSs6EBJ+DYQfbYhUEq3ID0HEEEI+88fZ4OGue885995z9t5n389rrb2yz97PPt9nk/C9z3328zxbtomIiGqNVF2BiIhIMo6IqIUk44iIGkgyjoiogSTjiIgaSDKOiKiBJOOIiClIWiXpfknbJZ0zyfn5kq4ozm+WtKQ4vkTSTyTdUWyf6hZrXg+VeQ2wGlhUHNoFjNu+dzo3FRExTCSNApcAJwA7gS2Sxm3f01ZsLfC47SMkrQEuAN5dnHvA9tG9xuvYMpb0e8AGQMAtxSbgC5P9lIiIaJAVwHbbO2zvoZULV08osxq4rNjfCKyUpJkE69YyXgscZfvZ9oOSPg5sAz422UWSxoAxgE9++Oxjxk4+aSZ1m7m9e8qN97yR0WriPlvR/c4/sJq4u39cfsz9Dyg/JsC8rr+8DsaeZyoJO/KGfzWjRNbu/Xppz9OKP82PzqTIVYV1ttcV+4uAh9rO7QSOnfAVL5SxvVfSE8Chxbmlkm4HngT+wPaNnerS7W96H/AK4B8mHD+sODep4mbWAey7/a8z3zoiaqk9V/XZ94BX2n5U0jHAX0o6yvaTU13QLRmfDfyNpL/npz8hXgkcAZzVjxpHRPRTH0cl7AIOb/u8uDg2WZmdkuYBC4FH3Vr05xkA27dKegB4FbB1qmAdk7HtayW9ilbfSfsDvC22n+v5liIiSjJvZl22k9kCHClpKa28twY4fUKZceAM4CbgFOAG25b0MuAx289J+ufAkcCOjvXuVhvb+4Cbp30bEREVGOlTLi76gM8CrgNGgUttb5N0PrDV9jiwHrhc0nbgMVoJG+B44HxJz9Lq0n2/7cc6xavo6UBExGD0c/KE7U3ApgnHzm3b3w2cOsl1VwJXTidWknFENMpI/7opSpVkHBGNMqzTipOMI6JR+tVnXLYk44holNF0U0REVC/dFBERNZBuioiIGkjLeCpVLdozb//SQ170xnd3LzQAv3vj5yuJW9nf7Wj5CzJd+7bfKD0mwKpr1lcS10/9sJK4/ZChbXVSQSKOiHqYN5y5uKHJOCLmrHRTRETUwAjD2TROMo6IRsloioiIGkg3RUREDaRlHBFRA31cXL5UScYR0SjppoiIqIF0U0RE1ECGtkVE1EBaxhERNTCaZBwRUb1h7aaY8YNHSe/tZ0UiIvphRL1vdTKbUSAfneqEpDFJWyVtXXfVNbMIERExPSPT2OqkYzeFpP831Sng5VNdZ3sdsA5g35ZNnnHtIiKmqWYN3p516zN+OfA24PEJxwX834HUKCJiFpq6uPzVwEG275h4QtLXB1KjiIhZqFv3Q686JmPbazucO73/1YmImJ3hbBdnaFtENIwa2k0RETFUhjMVJxlHRMM0ss84ImLYDGkvRZJxRDTLsE6HTjKOiEYZzlScZBwRDVO3NSd6lWQcEY2iIW0bDz4Zj4wOPMREX3/bGaXHBDj7+k9XEvdvT/wPlcR989XV3O9XTnxf6TFXf+e+0mMC3PnzR1cS97VfWV9J3H7oZyqWtAr4BDAKfMb2xyacnw98DjgGeBR4t+0H286/ErgHOM/2hZ1iDesokIiISfVrCU1Jo8AlwInAMuA0ScsmFFsLPG77COAi4IIJ5z8O9LR0ZZJxRDTKCOp562IFsN32Dtt7gA3A6gllVgOXFfsbgZUqpgBKeifwbWBbb/WOiGgQTWdrW3u92MbavmoR8FDb553FMSYrY3sv8ARwqKSDgN+jw7rvE+UBXkQ0ynQmfbSvvd5n5wEX2X6q17UykowjolH6+ABvF3B42+fFxbHJyuyUNA9YSOtB3rHAKZL+G3AIsE/SbtsXTxUsyTgiGqWPQ9u2AEdKWkor6a4BJi4dPA6cAdwEnALcYNvAm16oj3Qe8FSnRAxJxhHRMKN9ysW290o6C7iO1tC2S21vk3Q+sNX2OLAeuFzSduAxWgl7RpKMI6JR+jnO2PYmYNOEY+e27e8GTu3yHef1EivJOCIaJTPwIiJqIEtoRkTUwLBOnuhab0mvkbSyGMTcfnzV4KoVETEz05n0UScdk7Gk3wH+D/DbwN2S2qcC/ukgKxYRMRMjUs9bnXRrGf8mcIztdwJvAf5Q0geLc1PeSfsUw3Vf3jRVsYiIvhvWlnG3PuMR208B2H5Q0luAjZJ+jg730j7FcN+t17lPdY2I6KrX6cd1061l/H1JLyyoWiTmtwM/C7x2kBWLiJiJfi2hWbZuLeP3AHvbDxQrE71HUjUri0dEdKC6ZdkedUzGtnd2OPet/lcnImJ2RoZ0bFvGGUdEowxrn3GScUQ0ypDm4iTjiGiWtIwjImpgSHNxknFENEvdZtb1Ksk4IhplpIlD2yIiho0ytC0ionrD+gBPrXfnDc6V/+Sflb42xb+5+xtlh2zZ74Bq4lZkyy8eV0ncN9z61dJjauHLS48JsG/HHZXE3fCmf1tJ3NN/+PCsM+kDy47sOef8i3v+vjaZOy3jiGiUYW0ZJxlHRKMMaS5OMo6IZhnNaIqIiOqlmyIiogaGNBcnGUdEsyQZR0TUQCMXl4+IGDZ5gBcRUQPppoiIqIGMpoiIqIEhzcXdk7GkFYBtb5G0DFgF3Gd708BrFxExTY1sGUv6CHAiME/SV4Fjga8B50h6ne0/KaGOERE9G9Jc3LVlfApwNDAf+Edgse0nJV0IbAYmTcaSxoAxgDMXHMwJ8xf0r8YRER2MjA5nNu62DPNe28/Zfhp4wPaTALZ/Auyb6iLb62wvt708iTgiyiSp561OuiXjPZKez6bHPH9Q0kI6JOOIiMqMqPetC0mrJN0vabukcyY5P1/SFcX5zZKWFMdXSLqj2O6U9K6u1e5y/viiVYzt9uS7H3BG1zuJiCib1PvW8Ws0ClxC67nZMuC0YhBDu7XA47aPAC4CLiiO3w0st300rUEPn5bUsVu4YzK2/cwUxx+xfVfHO4mIqEAfuylWANtt77C9B9gArJ5QZjVwWbG/EVgpSbaftr23OH4A0PXtI0P66r6IiCmMjvS8SRqTtLVtG2v7pkXAQ22fdxbHmKxMkXyfAA4FkHSspG3AXcD725LzpDLpIyIaZToLBdleB6wbRD1sbwaOkvTzwGWSrrG9e6ryaRlHRLP0qc8Y2AUc3vZ5cXFs0jJFn/BC4NH2ArbvBZ4CfqFTsCTjiGgUjajnrYstwJGSlkraH1gDjE8oM85PBzOcAtxg28U18wAk/RzwGuDBTsHSTRERzdKn8cO290o6C7gOGAUutb1N0vnAVtvjwHrgcknbgcdoJWyA42jNVH6W1jDg37L9SKd4ScYR0Sx9XM+4WINn04Rj57bt7wZOneS6y4HLpxMryTgiGkWjw9n7mmQcEc1Ss2nOvZLddSzyrOzbsmmwASbxjZPGuhcagDd99tzuhZpk/gGVhPX4xtJjzvuvl3UvNAB/s6TjA/iB+dW//lwlcUd+aeWsM+nTJ7+x55yz4Mpv1SZzp2UcEc0ypC3jJOOIaJS8HToiog7SMo6IqF5GU0RE1EG6KSIiaiDdFBER1avb65R6lWQcEc2SboqIiOrlAV5ERB2kmyIionrDOulj2u15SdVMWo+I6EX/3vRRqo4tY0kTV7UX8KuSDgGw/Y5BVSwiYkYa2jJeDDwJfBz478X2o7b9SbW/cXXdVdf0q64REV1J6nmrk259xsuBDwIfBv6T7Tsk/cT233a6qP2Nq1UsoRkRc1gTR1PY3gdcJOlLxZ/f73ZNRESlatbi7VVPidX2TuBUSSfR6raIiKinJifj59n+K+CvBlSXiIjZG2lgN0VExNCZCy3jiIjaSzKOiKiB0dGqazAjScYR0SxpGUdE1ECScUREDSQZR0TUQIa2RUTUQJLxFLxv4CEmOv7qT5UeE+Drbz+zkrhvOvmXKok78uZfqySujntz6TH/fNGrS48JcPqNX6wkLrufriZuP6SbIiKiekrLOCKiBtIyjoiogSTjiIgaSDKOiKiBIZ0OPZw93RERU+njC0klrZJ0v6Ttks6Z5Px8SVcU5zdLWlIcP0HSrZLuKv7sOvQoyTgimqVPyVjSKHAJcCKwDDhN0rIJxdYCj9s+ArgIuKA4/gjwr22/FjgDuLxbtZOMI6JZRkZ63zpbAWy3vcP2HmADsHpCmdXAZcX+RmClJNm+3fZ3i+PbgAMlze9Y7WndZERE3U2jZdz+JvtiG2v7pkXAQ22fdxbHmKyM7b3AE8ChE8qcDNxm+5lO1c4DvIholmmMpmh/k/1gqqKjaHVdvLVb2STjiGiW/o2m2AUc3vZ5cXFssjI7Jc0DFgKPAkhaDFwFvMf2A92CTaubQtJxkj4kqWuWj4ioRP9GU2wBjpS0VNL+wBpgfEKZcVoP6ABOAW6wbUmH0Hp58zm2v9VLtTsmY0m3tO3/JnAxcDDwkcmGeUREVK5PybjoAz4LuA64F/ii7W2Szpf0jqLYeuBQSduBDwHP58WzgCOAcyXdUWz/tFO8bt0U+7XtjwEn2P6BpAuBm4GPTXZR0Qk+BvDJcz7A2LtWdQkTEdEnfVwoyPYmYNOEY+e27e8GTp3kuj8G/ng6sbol4xFJP0OrBS3bPygC/VjS3qkuau8U33fL1Z5OhSIiZqWh06EXArcCAizpMNvfk3RQcSwiol5GhnM6dMdkbHvJFKf2Ae/qe20iImZrZDjbiTMa2mb7aeDbfa5LRMTsaTjnsmWccUQ0S0P7jCMihkteuxQRUQNpGUdE1EATR1NERAyddFNERNRAuikiImogQ9siImpgLk36mJZ9FSxNsX/Ht5sMzFuu/WwlcW896X2VxH3Dn66vJO6WZceWHvP0r/5Z6TEB2P10NXGH9CEYMLR1T8s4Ipol3RQRETWQboqIiBrIaIqIiBpIN0VERA2kmyIiogYymiIiogbSTRERUQPppoiIqIG0jCMiaiBD2yIiaiBLaEZE1MCQjqbo+CNE0rGSXlrsHyjpo5K+IukCSQvLqWJExDRIvW810q09fynw/LJRnwAWAhcUx6ZcokzSmKStkrau+8tr+1LRiIiejIz0vtVIt26KEdt7i/3ltl9f7H9T0h1TXWR7HbAOYN/NX6lgDc2ImLNq1uLtVbcfDXdLem+xf6ek5QCSXgU8O9CaRUTMhEZ632qkW8v4fcAnJP0B8Ahwk6SHgIeKcxER9TKkD/A6JmPbTwC/UTzEW1qU32n7+2VULiJi2po8A8/2k8CdA65LRMTs1az7oVcZZxwRzTKkD/CSjCOiWYa0ZTyctY6ImIKknrcevmuVpPslbZd0ziTn50u6oji/WdKS4vihkr4m6SlJF/dS7yTjiGiWkXm9bx1IGgUuAU4ElgGnSVo2odha4HHbRwAX0ZoUB7Ab+EPgP/Zc7V4LRkQMhRH1vnW2Athue4ftPcAGYPWEMquBy4r9jcBKSbL9Y9vfpJWUe6t2rwUjIobCNCZ9tC/dUGxjbd+0iNaciuftLI4xWZlitvITwKEzqXYe4EVEs0xjNEX70g1VSzKOiGbp32iKXcDhbZ8XF8cmK7NT0jxai6k9OpNgg0/G8yrI93v3lB8TKhtSc8zV/7uSuLcsO7aSuMvHP1V+0OeeKz8mwGhFU3ur+P+2X/o3zngLcKSkpbSS7hrg9AllxoEzgJuAU4AbbM9ocbQh/i8eETGJPv0As71X0lnAdcAocKntbZLOB7baHgfWA5dL2g48RithAyDpQeClwP6S3gm81fY9U8VLMo6IZunjb6i2NwGbJhw7t21/N3DqFNcumU6sJOOIaJZMh46IqIEhnQ6dZBwRzZKWcUREDYwOZ1obzlpHREyhlwWA6ijJOCKaJX3GERE1kJZxREQNpGUcEVEDQ9oy7vgjRNLvSDq8U5mIiFoZHe19q5Fu7fk/AjZLulHSb0l6WRmVioiYsWmsZ1wn3Wqzg9aycX8EHAPcI+laSWdIOniqi9oXbF735Wv6WN2IiC6k3rca6dZnbNv7gOuB6yXtR+t9UKcBFwKTtpTbF2zet/WaGS0nFxExM/VKsr3qloxfdFe2n6W1fue4pAUDq1VExEzVrMXbq27J+N1TnbD9dJ/rEhExe01Mxrb/rqyKRET0Rc0ezPUq44wjolmGs2GcZBwRTTOc2TjJOCKapYl9xhERQyfJOCKiBvIALyKiDtIyjoioXropIiJqIMl4Cq5gaYqRei2N11THfP5j1QR+7rnyY+4/v/yYAHueqSbu3r3VxO2LJOOIiMrlhaQREXWQ0RQRETWQlnFERA0kGUdE1EGScURE9dIyjoiogeHMxUnGEdEwGU0REVED6aaIiKiD4UzGw9mej4iYitT71vWrtErS/ZK2SzpnkvPzJV1RnN8saUnbud8vjt8v6W3dYnVMxpL2l/QeSb9efD5d0sWSPiBpv653EhFRtj4lY0mjwCXAicAy4DRJyyYUWws8bvsI4CLgguLaZcAa4ChgFfC/iu+bUrduis8WZRZIOgM4CPgysBJYAZzR5fqIiHL17wHeCmC77R0AkjYAq4F72sqsBs4r9jcCF6u1OMZqYIPtZ4BvS9pefN9NUwXrloxfa/sXJc0DdgGvsP2cpD8H7pzqIkljwFjx8Uzb67rEmfJ7ZnrtbMyluHPpXuda3Ll0ry+yYGHPncYTchXAura6LwIeaju3Ezh2wle8UMb2XklPAIcWx2+ecO2iTnXp9iNkRNL+wMHAAmBhcXw+MGU3he11tpcX22z+Usa6FxmIuRR3Lt3rXIs7l+51Ribkqtnmq1np1jJeD9wHjAIfBr4kaQfwy8CGAdctIqJKu4DD2z4vLo5NVmZn0YOwEHi0x2tfpGPL2PZFwHHAr9j+n8DJwHXAWtsf7XorERHDawtwpKSlRQ/BGmB8Qplxfvrs7BTgBtsujq8pRlssBY4EbukUrOs4Y9vfbdv/Ia1O6rJU9SvDXIo7l+51rsWdS/fad0Uf8Fm0GqCjwKW2t0k6H9hqe5xW78HlxQO6x2glbIpyX6T1sG8v8AHbHV9RI1fxWqSIiHiRTPqIiKiBJOOIiBqobTLuNg1xQDEvlfSwpLvLiFfEPFzS1yTdI2mbpA+WFPcASbdIurOIW9oDWUmjkm6XdHVZMYu4D0q6S9IdkraWFPMQSRsl3SfpXkm/UkLMVxf3+Pz2pKSzBx23iP27xb+nuyV9QdIBZcRtglr2GRfTBv8OOIHWYOktwGm27+l44ezjHg88BXzO9i8MMlZbzMOAw2zfJulg4FbgnSXcq4CX2H6qmNr+TeCDtm/ucmk/Yn8IWA681PbbBx2vLe6DwHLbj5QY8zLgRtufKZ7ILygehJcVf5TWkKpjbf/DgGMtovXvaJntnxQPsDbZ/rNBxm2KuraMX5iGaHsPrTHNqwcd1PY3aD0RLY3t79m+rdj/EXAvXWbq9CmubT9VfNyv2Ab+k1nSYuAk4DODjlU1SQuB42k9ccf2njITcWEl8MCgE3GbecCBxZjbBcB3u5SPQl2T8WTTEAeeoKpWrPj0OmBzSfFGJd0BPAx81XYZcf8H8J+BfSXEmsjA9ZJuLabBDtpS4AfAZ4tumc9IekkJcdutAb5QRiDbu4ALge8A3wOesH19GbGboK7JeM6RdBBwJXC27SfLiGn7OdtH05odtELSQLtmJL0deNj2rYOM08Fxtl9PaxWuDxTdUoM0D3g98EnbrwN+DJTy/ANaqy4C7wC+VFK8n6H1G+xS4BXASyT9uzJiN0Fdk/G0pxIOs6LP9krgL2x/uez4xa/OX6O11N8gvRF4R9F3uwH4tWLRqVIULTdsPwxcRas7bJB2AjvbfuPYSCs5l+VE4Dbb3y8p3q8D37b9A9vP0lrh8V+WFHvo1TUZ9zINsRGKB2nrgXttf7zEuC+TdEixfyCth6X3DTKm7d+3vdj2Elp/pzfYLqXlJOklxQNSiq6CtwIDHTVj+x+BhyS9uji0khcvvzhop1FSF0XhO8AvS1pQ/LteSesZSPSglq9dmmoa4qDjSvoC8BbgZyXtBD5ie/2Aw74R+PfAXUX/LcB/sb1pwHEPAy4rnraPAF+0XepQs5K9HLiqlSOYB3ze9rUlxP1t4C+KRsUO4L0lxHz+B84JwJllxAOwvVnSRuA2WlOAb6chU6PLUMuhbRERc01duykiIuaUJOOIiBpIMo6IqIEk44iIGkgyjoiogSTjiIgaSDKOiKiB/w/qUo35oXt00QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"Predicted data uncertainty:\")\n",
    "sns.heatmap(var_preds.reshape([num_cat, num_cat]), cmap=\"Reds\", vmin=0, vmax=0.05)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, we know how to estimate noise in the data. But how to measure knowledge uncertainty coming from the lack of training data in a particular region? What to do if we want to detect outliers?  Estimating knowledge uncertainty requires an ensemble of models. If all the models give similar predictions, then there is low knowledge uncertainty, and if they give diverse predictions (strongly disagree with each other), then there is high knowledge uncertainty.\n",
    "Let’s generate an ensemble:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ensemble(train_pool, val_pool, num_samples=10, iters=1000, lr=0.2):\n",
    "    ens_preds = []\n",
    "    for seed in range(num_samples):\n",
    "        model = CatBoostRegressor(random_seed=seed, iterations=iters, learning_rate=lr, \n",
    "                                  loss_function='RMSEWithUncertainty', posterior_sampling=True,\n",
    "                                  verbose=False)\n",
    "        model.fit(train_pool, eval_set=val_pool)\n",
    "        ens_preds.append(model.predict(test))\n",
    "    return ens_preds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the models are generated using the option posterior_sampling, since this allows the obtained (random) predictions to be nicely distributed (with good theoretical properties).\n",
    "\n",
    "For regression, knowledge uncertainty can be obtained by measuring the variance of the mean across multiple models. Below we also plot the average estimated data uncertainty (average predicted variance for the modes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Knowledge uncertainty via ensemble:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD8CAYAAABekO4JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAZEklEQVR4nO3dfZBddZ3n8fenOwQJD3EWHwoTxqSWOApMlQoV3MGiHGMwzDCGKaGIwwo6jL2WZtSd2oe4M2CJ7tZStTs+jJRrNJGHVUCDlL1DBGXBdWXWkABxIAHcJoJ0YHmeIESMgc/+cU/kptN9z+30vbfPPffzqjrV556n7+9UOt/+3e/5nXNkm4iIqK6h2W5ARES0lkQdEVFxSdQRERWXRB0RUXFJ1BERFZdEHRFRcUnUEREVN6dsA0lvBFYCC4pFO4FR2/d2s2EREdHQskct6d8D1wACbi8mAVdLWtP95kVEhFrdmSjpZ8AJtn8zYflcYJvtJVPsNwKMAHzl7z5/0siff6BjDY6IGps3XzM9xId1VNu3W/83PzvjeL1QVvp4CXgd8NCE5ccU6yZley2wFoDdu3KPekTEDJQl6k8A/1PS/wUeLpb9LnAcsLqbDYuIOBh1HCHRMlHbvlHSG4Cl7H8xcbPtF7vduIiI6ZqjvqhmTEvpqA/bLwE/6UFbIiJmbKh+ebo8UUdE9JOBK31ERPSboUEsfURE9JP0qCMiKi416oiIihtO6SMiotpS+oiIqLiUPiIiKi496oiIisvwvIiIiptTvzydRB0R9VLH0kcdzykiBtgQansqI2mFpPsljU32shRJh0q6tli/SdKiYvlSSVuL6aeS/rTdY05+ThERNTKk9qdWJA0DlwFnAMcD75N0/ITNLgSesX0c8Dng0mL5PcDJtt8MrAC+ImlOm8c88JzaPfmIiH4wNI2pxFJgzPYO23tovJZw5YRtVgJXFPMbgGWSZHu37b3F8lcA+16g0s4xJz2niIjamE6PWtKIpC1N00jToRbw8gtTAMZ5+bn8B2xTJOZdwNEAkk6RtA24G/hwsb6dYx4gFxMjolam8+KA/V4b2GG2NwEnSHoTcIWk7x3ssdKjjoha6WDpYydwbNPnhcWySbeRNAeYDzzVvIHte4HngBPbPOak5xQRURudupgIbAaWSFosaS6wChidsM0ocEExfzZwi20X+8wBkPR64I3Ag20e8wApfURErbQz7K4dtvdKWg3cBAwD621vk3QJsMX2KLAOuErSGPA0jcQL8HZgjaTfAC8BH7H9JMBkxyxri2yXbTMzu3d1OUBE1Ma8+TPOsuuOenXbOefCZ5/oi/sY06OOiFoZ7ovUOz1J1BFRK50qfVTJQV9MlPTBTjYkIqITOngxsTJmMurj01OtaB5Evnb95TMIERExPR0cnlcZLUsfkv5xqlXAa6fab79B5LmYGBE91Ecd5baV1ahfC7wbeGbCcgH/0JUWRUTMwCC+OODvgSNsb524QtIPu9KiiIgZ6KeSRrtaJmrbF7ZY92edb05ExMzUrz+d4XkRUTMawNJHRERfqV+aTqKOiJoZuBp1RES/qWHlI4k6IuqljreQJ1FHRK3UL00nUUdEzfTTMzzalUQdEbWiGvapk6g76MV/+O6sxB3+g9K3zccM7fno2bMSd+5lG2Ylbj+rX5pOoo6ImknpIyKi4jLqIyKi4uqXppOoI6JmcsNLRETF1TBPJ1FHRL1keF5ERMUN1y9P1/JBUxExwDSNqfRY0gpJ90sak7RmkvWHSrq2WL9J0qJi+XJJd0i6u/j5zqZ9flgcc2sxvaasHelRR0StdKr0IWkYuAxYDowDmyWN2t7etNmFwDO2j5O0CrgUOBd4EvgT249IOhG4CVjQtN95tre025b0qCOiVqT2pxJLgTHbO2zvAa4BJt4GvBK4opjfACyTJNt32X6kWL4NOEzSoQd7TknUEVErQ9OYJI1I2tI0jTQdagHwcNPncfbvFe+3je29wC7g6AnbvBe40/avm5Z9vSh7XKQ23h1WWvqQ9MaiMZtsP9e0fIXtG8v2j4jopekUPmyvBdZ2rS3SCTTKIac3LT7P9k5JRwLXAe8Hrmx1nJY9akkfA74L/CVwj6Tmbv9/OpiGR0R005DU9lRiJ3Bs0+eFxbJJt5E0B5gPPFV8XghcD5xv+4F9O9jeWfz8JfBNGiWW1udUsv5DwEm2zwLeAVwk6ePFuinPsvnrxNr1l5e1ISKiYzo46mMzsETSYklzgVXA6IRtRoELivmzgVtsW9IrgRuANbZv+23bpDmSXlXMHwKcCdxT1pCy0sfQvnKH7QclvQPYIOn1tDjP/b5O7N7lskZERHRKGyXfttjeK2k1jREbw8B629skXQJssT0KrAOukjQGPE0jmQOsBo4DLpZ0cbHsdOB54KYiSQ8DNwNfLWtLWaJ+TNKbbW8tGv6cpDOB9cDvt3/KERG90cnHnNreCGycsOzipvkXgHMm2e+zwGenOOxJ021HWaI+H9g7oQF7gfMlfWW6wSIiuk01fCB1y0Rte7zFutumWhcRMVuGajjoOHcmRkStdKpGXSVJ1BFRKzXM00nUEVEv6VFHRFRcDfN0EnVE1Esbdxz2nSTqiKiVoUEbnhcR0W+U4XkREdWWi4l9Yu/FF85K3L/5ws2zEhc+NitR//Ouh2Yl7ocPP7Z8ow47858d3vOYAG/ZesqsxF1w26ZZidsJNczT9UzUETG40qOOiKi4GubpJOqIqJfhjPqIiKi2lD4iIiquhnk6iToi6iWJOiKi4gbuxQEREf0mFxMjIioupY+IiIrLqI+IiIqrYZ4uT9SSlgK2vVnS8cAK4L7iNeoREZUycD1qSZ8CzgDmSPoBcApwK7BG0lts/8cetDEiom01zNOUPbn1bOBU4DTgo8BZtj8DvBs4d6qdJI1I2iJpy9r1l3eqrRERpYaG1fZURtIKSfdLGpO0ZpL1h0q6tli/SdKiYvlySXdIurv4+c6mfU4qlo9J+qLa+ApQVvrYa/tFYLekB2w/C2D7V5Jemmon22uBtQDs3uWyRkREdEqnSh+ShoHLgOXAOLBZ0qjt7U2bXQg8Y/s4SauAS2l0Yp8E/sT2I5JOBG4CFhT7fBn4ELAJ2EijnPy9Vm0p61HvkTSvmD+p6QTmA1Mm6oiIWTOk9qfWlgJjtnfY3gNcA6ycsM1K4IpifgOwTJJs32X7kWL5NuCwovd9DHCU7Z/YNnAlcFbpKZWsP832bgDbzYn5EOCCsoNHRPSc1PbUXKYtppGmIy0AHm76PM7LveIDtrG9F9gFHD1hm/cCd9r+dbH9eMkxD9Cy9FEceLLlT9Lo2kdEVMp0Sh/7lWm705YTaJRDTp/JcTKOOiLqZbhjb7fdCTS/921hsWyybcYlzQHmA08BSFoIXA+cb/uBpu0XlhzzADV8X29EDDINqe2pxGZgiaTFkuYCq4DRCduM8nIZ+GzgFtuW9ErgBmCN7dv2bWz7UeBZSW8rRnucD3y3rCFJ1BFRL9OoUbdS1JxX0xixcS/wLdvbJF0i6T3FZuuAoyWNAX8F7BvCtxo4DrhY0tZiek2x7iPA14Ax4AFKRnxASh8RUTOdfMxpcQf2xgnLLm6afwE4Z5L9Pgt8dopjbgFOnE47kqgjol5qeGtiEnVE1EueRx0RUW3q3KiPykiijoh6qWHpQ427GLtogJ71seuPTpvtJvTUEX906uwEXvKmnod87ovreh4TYP4N/2tW4s6aefNnnGV3v/fUtnPOvOtu64usnh51RNRLDXvUSdQRUSt5C3lERNWlRx0RUW0Z9RERUXUpfUREVFxKHxER1TZwbyGPiOg7KX1ERFRbLiZGRFRdSh8REdVWxxtepv0dQdKV3WhIRERHdOgNL1XSskctaeL7wQT8YfE+MGy/58C9IiJmUQ171GWlj4XAdhrv9zKNRH0y8F9b7SRpBBgB+MrffZ6RP//AjBsaEdGOQRyedzLwceCvgX9re6ukX9lu+exF22uBtcBAPeY0Iipg0EZ92H4J+Jykbxc/HyvbJyJiVg1gjxoA2+PAOZL+GHi2u02KiJiBQU3U+9i+AbihS22JiJi5oQErfURE9J0a9qjr96cnIgZbB8dRS1oh6X5JY5LWTLL+UEnXFus3SVpULD9a0q2SnpP0pQn7/LA45tZiek1ZO9Kjjoh6GR7uyGEkDQOXAcuBcWCzpFHb25s2uxB4xvZxklYBlwLnAi8AFwEnFtNE59ne0m5b0qOOiHrpXI96KTBme4ftPcA1wMoJ26wErijmNwDLJMn287Z/TCNhz1gSdUTUyzQStaQRSVuappGmIy0AHm76PF4sY7JtbO8FdgFHt9HKrxdlj4vUxh06KX1ERL1M42Lifjfn9c55tndKOhK4Dng/0PIZSulRR0S9DA21P7W2Ezi26fPCYtmk20iaA8wHnmp1UNs7i5+/BL5Jo8TS+pTKNoiI6CudS9SbgSWSFkuaC6wCJj6obhS4oJg/G7jF9pSPzZA0R9KrivlDgDOBe8oaktJHB83f+KNZifuzN580K3GPWPTPZyUud9/V85Dzb2j5eJuokg6No7a9V9Jq4CZgGFhve5ukS4AttkeBdcBVksaAp2kk86IZehA4Cpgr6SzgdOAh4KYiSQ8DNwNfLWtLEnVE1Io6eGei7Y3AxgnLLm6afwE4Z4p9F01x2Gn3rJKoI6JeanhnYhJ1RNRLEnVERMUlUUdEVFyHbiGvkiTqiKiX9KgjIiouiToiouLy4oCIiIpLjzoiouKSqCMiKm7QR31IejuNJz3dY/v73WlSRMQM1LBH3bLqLun2pvkPAV8CjgQ+Ndn7wyIiZl0H35lYFWWXRw9pmh8Bltv+NI2nQJ031U7Nb01Yu/7ymbcyIqJdnXvMaWWUlT6GJP0OjYQu208A2H5e0t6pdtrvrQm7d035bNaIiI7ro55yu8oS9XzgDkCAJR1j+1FJRxTLIiKqZWjALia2eJ7qS8Cfdrw1EREzNVS/PuRBDc+zvRv4eYfbEhExc+qf2nO7Mo46IuplAGvUERH9pY9Gc7QriToi6iU96oiIihu0UR8REX0npY+IiIpL6SMiouJqODyvfmcUEYNtSO1PJSStkHS/pLHJHkQn6VBJ1xbrN0laVCw/WtKtkp6T9KUJ+5wk6e5iny9K5V8B0qOugTdsvWNW4r64/jOzEnf433xuVuJGn+jQxURJw8BlwHJgHNgsadT29qbNLgSesX2cpFXApcC5wAvARcCJxdTsy8CHgE3ARmAF8L1WbUmPOiLqRUPtT60tBcZs77C9B7gGWDlhm5XAFcX8BmCZJNl+3vaPaSTsl5smHQMcZfsntg1cCZxV1pAk6oiol86VPhYADzd9Hi+WTbqN7b3ALuDokmOOlxzzwFMq2yAioq9M48UBzc/OL6aR2W7+ZFKjjoh6mcaoj/2enX+gncCxTZ8XFssm22Zc0hwaj4Z+qkXIncVxWh3zAOlRR0S9dK70sRlYImmxpLnAKmB0wjajwAXF/NnALUXteVK2HwWelfS2YrTH+cB3yxqSHnVE1EuHRn3Y3itpNXATMAyst71N0iXAFtujwDrgKkljwNM0kjkAkh4EjgLmSjoLOL0YMfIR4HLgMBqjPVqO+IAk6oiomw7e8GJ7I40hdM3LLm6afwE4Z4p9F02xfAsHDtlrKYk6Iuolb3iJiKi4Gt5CnkQdEfWShzJFRFRcHnMaEVFxNXxxQMs/PZJOkXRUMX+YpE9L+h+SLpU0vzdNjIiYhmncmdgvyr4jrAd2F/NfoHHXzaXFsq9PtVPzbZlr11/eiXZGRLRnaKj9qU+UlT6GigeNAJxs+63F/I8lbZ1qp/1uy9y9a8q7dCIiOq6PesrtKvuTco+kDxbzP5V0MoCkNwC/6WrLIiIORucec1oZZT3qvwC+IOlvgCeB/yPpYRqP9fuLbjcuImLaangxsWWitr0L+EBxQXFxsf247cd60biIiGkb1DsTbT8L/LTLbYmImLk+Kmm0K+OoI6JeangxMYk6IuolPeqIiGpTetQRERU3VL+0Vr8ziojBNqijPiIi+kZq1BERFZcadURExaVH3Semflt7twPPTthZ+sUcvuCTsxJ3Vrz04uzEreHt0F2XHnVERMUN1++PWxJ1RNRLSh8RERWX0kdERMXVsEddvzOKiMHWwXcmSloh6X5JY5LWTLL+UEnXFus3SVrUtO6TxfL7Jb27afmDku6WtFXSlnZOKT3qiKiX4c6kNUnDwGXAcmAc2Cxp1Pb2ps0uBJ6xfZykVTTeKXuupOOBVcAJwOuAmyW9wfa+4UN/aPvJdtuSHnVE1IqktqcSS4Ex2zts7wGuAVZO2GYlcEUxvwFYpsaBVwLX2P617Z8DY8XxDkoSdUTUyzTemShpRNKWpmmk6UgLaLx2cJ/xYhmTbVO8CHwXcHTJvga+L+mOCfGmlNJHRNTLNEZ92F4LrO1eYyb1dts7Jb0G+IGk+2z/qNUO6VFHRL107i3kO4Fjmz4vLJZNuo2kOcB84KlW+9re9/Nx4HraKIkkUUdEvXRu1MdmYImkxZLm0rg4ODphm1HggmL+bOAW2y6WrypGhSwGlgC3Szpc0pGNZupw4HTgnrKGtCx9SPoYcL3th1ttFxFRGR26hdz2XkmrgZuAYWC97W2SLgG22B4F1gFXSRoDnqaRzCm2+xawHdgLfNT2i5JeC1xfXMicA3zT9o1lbZFbPMBI0i7geeAB4Grg27afmNbZ7t7V+ycV5aFMvfHi3tmJ26HhV9OShzL1xrz5M76t0A/d3fZ/RL3+9/viNsay/+E7aNRWPgOcBGyXdKOkC/Z13yfTfCV17frLO9faiIgyHbzhpSrKetR32n5r0+dDgDOA9wHvsv3q0gjpUXdfetTdlx51b3SiR/2Lbe33qH/3hL7I1mW/8fudhO3f0CiSj0qa17VWRUQcrD7qKberLFGfO9UK27s73JaIiJkbtERt+2e9akhEREfU8Ol5uTMxIuqlfh3qJOqIqJv6Zeok6oiol0GrUUdE9J0k6oiIisvFxIiIqkuPOiKi2lL6iIiouCTqPjFr/1D1+wVpaTaeuTFbBu2ZG32tfv8PB+h/WkQMgjZeWtt3kqgjol4y6iMiouLSo46IqLgk6oiIqkuijoiotvSoIyIqrn55Ook6Imomoz4iIioupY+IiKpLoo6IqLZB61FLmgusAh6xfbOkPwP+ALgXWGv7Nz1oY0RE+2qYqGV76pXSN2gk83nAPwFHAN8BlhX7XlAaYfeuqQNERDSbN3/mWXY6OacT8XrB9pQT8I/FzznAY8Bw8Vn71k2x3wiwpZhGWsUoiX/Q+85kGqS4g3SugxZ3kM617lNZj/oe4K3A4cAvgNfbflrSK4C7bL+pM38upoy/xfbJ3Ywx6HEH6VwHLe4gnWvdlV1MXAfcBwwDfw18W9IO4G3ANV1uW0REUJKobX9O0rXF/COSrgTeBXzV9u29aGBExKArHZ5n+5Gm+X8CNnS1Rftb28NYgxp3kM510OIO0rnWWssadUREzL763RQfEVEzSdQRERVX2UQtaYWk+yWNSVrTo5jrJT1eDEvsCUnHSrpV0nZJ2yR9vEdxXyHpdkk/LeJ+uhdxi9jDku6S9Pe9ilnEfVDS3ZK2StrSo5ivlLRB0n2S7pX0L3oQ8/eKc9w3PSvpE92OW8T+18Xv0z2Sri6G8sYMVbJGLWkY+BmwHBgHNgPvs729y3FPA54DrrR9YjdjNcU8BjjG9p2SjgTuAM7qwbkKONz2c5IOAX4MfNz2T7oZt4j9V8DJwFG2z+x2vKa4DwIn236yhzGvAP637a8Vj2SYV1yU71X8YWAncIrth7ocawGN36Pjbf9K0reAjbYv72bcQVDVHvVSYMz2Dtt7aIzZXtntoLZ/BDzd7TgTYj5q+85i/pc0nqOyoAdxbfu54uMhxdT1v9qSFgJ/DHyt27Fmm6T5wGk07kfA9p5eJunCMuCBbifpJnOAwyTte/TEIyXbRxuqmqgXAA83fR6nB8lrtklaBLwF2NSjeMOStgKPAz+w3Yu4nwf+HfBSD2JNZOD7ku6QNNKDeIuBJ4CvF6Wer0k6vAdxm60Cru5FINs7gf9C4y7mR4Fdtr/fi9h1V9VEPXAkHQFcB3zC9rO9iGn7RdtvBhYCSyV1tdwj6Uzgcdt3dDNOC2+3/VbgDOCjRamrm+bQeATDl22/BXge6Mn1Fvjt0y/fA3y7R/F+h8Y338XA64DDJf3LXsSuu6om6p3AsU2fFxbLaqmoEV8HfMP2d3odv/g6fiuwosuhTgXeU9SKrwHeKem/dznmbxU9Pmw/DlxPo8TWTePAeNM3lQ00EnevnAHcafuxHsV7F/Bz20+48Qjk79B4LHLMUFUT9WZgiaTFTc/EHp3lNnVFcVFvHXCv7b/tYdxXS3plMX8YjQu393Uzpu1P2l5oexGNf9NbbPekxyXp8OJiLUX54XSgq6N7bP8/4GFJv1csWgZ09SLxBO+jR2WPwi+At0maV/xeL6NxzSVmqJJveLG9V9Jq4CYaD4Rab3tbt+NKuhp4B/AqSePAp2yv63LYU4H3A3cX9WKA/2B7Y5fjHgNcUYwKGAK+Zbunw+V67LXA9Y38wRzgm7Zv7EHcvwS+UXQ4dgAf7EHMfX+MlgP/qhfxAGxvkrQBuBPYC9xFbifviEoOz4uIiJdVtfQRERGFJOqIiIpLoo6IqLgk6oiIikuijoiouCTqiIiKS6KOiKi4/w+IeR4R/Txv2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average predicted data uncertainty:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAD8CAYAAACihcXDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAX9UlEQVR4nO3de5RdZXnH8e/vnJBAuAQbLQsSNOkCtUFbBBq8IIsawSCXYA2LwFKiTR0ooFiXrVgUAW1ruijULigwJdiQIkECtAOGWw22QiEk3ArhYoeAMhG5mxggJGGe/nF29HCcc5mZc/beZ8/vs9Ze7LMv53kPmfXMO89+3/coIjAzs2yVsm6AmZk5GZuZ5YKTsZlZDjgZm5nlgJOxmVkOOBmbmeWAk7GZWR2SZkt6XFK/pDOGOD9B0tXJ+ZWSpiXHp0l6TdIDyXZJs1jjWmjMu4E5wJTk0DqgLyIeHc6HMjPrJpLKwEXAocAAsEpSX0Q8UnXZAuDliNhL0jxgIXBccu6JiNi31XgNe8aSvgIsBQTck2wCrhrqt4SZWYHMBPojYm1EbKaSC+fUXDMHWJzsLwNmSdJIgjXrGS8A9omILdUHJZ0PrAG+PdRNknqAHoCLzzx9/54/+fhI2jZym19PN9424ydkEzerz7vd+GzibtmcfsxxTf+I7IxSOZu4b2zNJGzpjz4+okRW7WTt0vK04kv51UkkuSrRGxG9yf4U4OmqcwPAgTVv8etrImKrpPXA5OTcdEn3AxuAr0XEjxu1pdlP2CCwB/DTmuO7J+eGlHyYXoDB+271fGszy6XqXNVmzwBvj4gXJe0P/LukfSJiQ70bmiXjLwI/lPR//OY3xNuBvYDT2tFiM7N2auOohHXAnlWvpybHhrpmQNI4YBLwYlQW/XkdICLulfQE8E5gdb1gDZNxRNws6Z1UaifVD/BWRcQbLX8kM7OUjBtZyXYoq4C9JU2nkvfmASfUXNMHzAfuAuYCKyIiJL0NeCki3pD0e8DewNqG7W7WmogYBO4e9scwM8tAqU25OKkBnwbcApSByyNijaRzgdUR0QcsApZI6gdeopKwAQ4GzpW0hUpJ9+SIeKlRvIyeSpiZdUY7J09ExHJgec2xs6r2NwHHDnHftcC1w4nlZGxmhVJqX5kiVU7GZlYo3Tqt2MnYzAqlXTXjtDkZm1mhlF2mMDPLnssUZmY54DKFmVkOuGdcTxYLjpTLUE7/98w1sz6dekyAY//zikziZrJgT0Zunv2nmcSdfduSTOLG669mErcdPLQtTzJIxGaWD+O6MxcXNBmb2ZjlMoWZWQ6U6M6usZOxmRWKR1OYmeWAyxRmZjngnrGZWQ60cXH5VDkZm1mhuExhZpYDLlOYmeWAh7aZmeWAe8ZmZjlQdjI2M8tet5YpRvzgUdJn29kQM7N2KKn1LU9GMwrknHonJPVIWi1pde/1N40ihJnZ8JSGseVJwzKFpP+tdwrYrd59EdEL9AIMrloeI26dmdkw5azD27JmNePdgI8BL9ccF/A/HWmRmdkoFHVx+RuBnSLigdoTkn7UkRaZmY1C3soPrWqYjCNiQYNzJ7S/OWZmo9Od/WIPbTOzglFByxRmZl2lO1Oxk7GZFUwha8ZmZt2mS6sUTsZmVizdOh3aydjMCqU7U7GTsZkVTN7WnGiVk7GZFYq6tG/c+WRcKnc8RK0VH8tmQbnbf/laJnGnfLzu3JyO+uANl2QS98bDP5d6zAOm7Zp6TIDbD5ufSdw/vnVxJnHboZ2pWNJs4DtAGbgsIr5dc34CcAWwP/AicFxEPFV1/u3AI8DZEXFeo1jdOgrEzGxI7VpCU1IZuAg4HJgBHC9pRs1lC4CXI2Iv4AJgYc3584GWlq50MjazQimhlrcmZgL9EbE2IjYDS4E5NdfMAbb9GbEMmKVkCqCkY4AngTWttdvMrEA0nK1q7fVk66l6qynA01WvB5JjDHVNRGwF1gOTJe0EfIUG677X8gM8MyuU4Uz6qF57vc3OBi6IiI2trpXhZGxmhdLGB3jrgD2rXk9Njg11zYCkccAkKg/yDgTmSvp7YFdgUNKmiLiwXjAnYzMrlDYObVsF7C1pOpWkOw+oXTq4D5gP3AXMBVZERAAf/nV7pLOBjY0SMTgZm1nBlNuUiyNiq6TTgFuoDG27PCLWSDoXWB0RfcAiYImkfuAlKgl7RJyMzaxQ2jnOOCKWA8trjp1Vtb8JOLbJe5zdSiwnYzMrFM/AMzPLAS+haWaWA906eaJpuyW9W9KsZBBz9fHZnWuWmdnIDGfSR540TMaSvgD8B/B54GFJ1VMB/7aTDTMzG4mS1PKWJ816xp8D9o+IY4BDgK9LOj05V/eTVE8x7L1ueb3LzMzarlt7xs1qxqWI2AgQEU9JOgRYJukdNPgs1VMMB++9JdrUVjOzplqdfpw3zXrGz0rad9uLJDEfCbwVeG8nG2ZmNhLtWkIzbc16xicCW6sPJCsTnSjp0o61ysxshJS3LNuihsk4IgYanLuz/c0xMxudUpeObfM4YzMrlG6tGTsZm1mhdGkudjI2s2Jxz9jMLAe6NBc7GZtZseRtZl2rnIzNrFBKRRzaZmbWbeShbWZm2fMDvDp+cNhnOh3itxx+3H6pxwQ4ZML4TOJm5fYjTsok7hFXnJN6TP3ObqnHBNjtrh9mEvf6WZ/KJO4nX/rFqN+jS3Oxe8ZmVizuGZuZ5UCX5mInYzMrlrJHU5iZZc9lCjOzHOjSXOxkbGbF4mRsZpYDhVxc3sys2/gBnplZDrhMYWaWAx5NYWaWA12ai5snY0kzgYiIVZJmALOBxyJiecdbZ2Y2TIXsGUv6BnA4ME7SbcCBwO3AGZLeFxF/k0Ibzcxa1qW5uGnPeC6wLzAB+AUwNSI2SDoPWAkMmYwl9QA9AKdM3IXZ209sX4vNzBoolbszGzdLxlsj4g3gVUlPRMQGgIh4TdJgvZsiohfoBbhh8u7RttaamTXRrWWKZmvib5a0rVu7/7aDkiYBdZOxmVlmSmp9a0LSbEmPS+qXdMYQ5ydIujo5v1LStOT4TEkPJNuDkj7RtNlNzh8cEa8CRER18t0OmN/0k5iZpU1qfWv4NioDF1F5bjYDOD4ZxFBtAfByROwFXAAsTI4/DBwQEftSGfRwqaSGlYiGyTgiXq9z/IWIeKjhJzEzy4CklrcmZgL9EbE2IjYDS4E5NdfMARYn+8uAWZIUEa9GxNbk+PZA03Jtl351n5lZHeVSy5ukHkmrq7aeqneaAjxd9XogOcZQ1yTJdz0wGUDSgZLWAA8BJ1cl5yF50oeZFcpwFgqqHmzQbhGxEthH0u8DiyXdFBGb6l3vnrGZFUubasbAOmDPqtdTk2NDXpPUhCcBL1ZfEBGPAhuB9zQK5mRsZoWiklremlgF7C1puqTxwDygr+aaPn4zmGEusCIiIrlnHICkdwDvBp5qFMxlCjMrljaNM46IrZJOA24BysDlEbFG0rnA6ojoAxYBSyT1Ay9RSdgAB1GZqbyFyjDgUyLihUbxnIzNrFjauJ5xsgbP8ppjZ1XtbwKOHeK+JcCS4cRyMjazQlG5O6uvTsZmVixdOh2648n4iJv+pdMhfsudR52cekyAK5/bkEncrHztD/fIJO6pR3459ZgXnnl06jEB7uj9USZxP/GjqzKJ2w7qzo6xe8ZmVjDuGZuZZc/fDm1mlgfuGZuZZc+jKczM8sBlCjOzHHCZwswse936tUtOxmZWLC5TmJllzw/wzMzywGUKM7Psdeukj2H35yVd0YmGmJm1Rfu+6SNVDXvGkmpXtRfwx5J2BYiIbFZPMTOrp6A946nABuB84B+S7VdV+0Oq/sbV3utvbldbzcyaktTylifNasYHAKcDZwJ/GREPSHotIv6r0U3V37g6eM+N0ZaWmpm1ooijKSJiELhA0jXJf59tdo+ZWaZy1uNtVUuJNSIGgGMlHUGlbGFmlk9FTsbbRMQPgB90qC1mZqNXKmCZwsys64yFnrGZWe45GZuZ5UC5nHULRsTJ2MyKxT1jM7MccDI2M8sBJ2Mzsxzw0DYzsxxwMq5jMP2lKT50wyWpxwQoHfXnmcRd8tz6TOLu/JbtM4m78MPTUo953UW3pR4TYO6KKzOJy6ZXsonbDi5TmJllT+4Zm5nlgHvGZmY54GRsZpYDTsZmZjnQpdOhu7PSbWZWTxu/kFTSbEmPS+qXdMYQ5ydIujo5v1LStOT4oZLulfRQ8t+PNIvlZGxmxdKmZCypDFwEHA7MAI6XNKPmsgXAyxGxF3ABsDA5/gJwVES8F5gPLGnWbCdjMyuWUqn1rbGZQH9ErI2IzcBSYE7NNXOAxcn+MmCWJEXE/RHx8+T4GmAHSRMaNntYH9LMLO+G0TOu/ib7ZOupeqcpwNNVrweSYwx1TURsBdYDk2uu+SRwX0S83qjZfoBnZsUyjNEU1d9k35mmaB8qpYvDml3rZGxmxdK+0RTrgD2rXk9Njg11zYCkccAk4EUASVOB64ETI+KJZsGGVaaQdJCkL0lqmuXNzDLRvtEUq4C9JU2XNB6YB/TVXNNH5QEdwFxgRUSEpF2pfHnzGRFxZyvNbpiMJd1Ttf854EJgZ+AbQw3zMDPLXJuScVIDPg24BXgU+H5ErJF0rqSjk8sWAZMl9QNfArblxdOAvYCzJD2QbL/bKF6zMsV2Vfs9wKER8byk84C7gW8PdVNSBO8BuPgrp9BzzOwmYczM2qSNCwVFxHJgec2xs6r2NwHHDnHft4BvDSdWs2RckvQWKj1oRcTzSaBXJG2td1N1UXzw7hvSX0PTzMaugk6HngTcCwgISbtHxDOSdkqOmZnlS6k7p0M3TMYRMa3OqUHgE21vjZnZaJW6s584oqFtEfEq8GSb22JmNnrqzrlsHmdsZsVS0JqxmVl38dcumZnlgHvGZmY5UMTRFGZmXcdlCjOzHHCZwswsBzy0zcwsB8bSpI9hyeJ/zHbbpx8T+MCNl2YSt3zUSZnE3WnOIZnEvfPvlqUec+6t3009JgCDb2QTt7xd82vyyg/wzMxywGUKM7MccJnCzCwHPJrCzCwHXKYwM8sBlynMzHLAoynMzHLAZQozsxxwmcLMLAfcMzYzywEPbTMzywEvoWlmlgNdOpqi4a8QSQdK2iXZ30HSOZJukLRQ0qR0mmhmNgxS61uONOvPXw68mux/B5gELEyO1V3GSlKPpNWSVvdef3NbGmpm1pJSqfUtR5qVKUoRsTXZPyAi9kv275D0QL2bIqIX6AUYvOfGGH0zzcxalLMeb6ua/Wp4WNJnk/0HJR0AIOmdwJaOtszMbCRUan3LkWY94z8DviPpa8ALwF2SngaeTs6ZmeVLlz7Aa5iMI2I98JnkId705PqBiHg2jcaZmQ1bkWfgRcQG4MEOt8XMbPRyVn5olccZm1mxdOkDPCdjMyuWLu0Zd2erzczqkNTy1sJ7zZb0uKR+SWcMcX6CpKuT8yslTUuOT5Z0u6SNki5spd1OxmZWLKVxrW8NSCoDFwGHAzOA4yXNqLlsAfByROwFXEBlUhzAJuDrwJdbbnarF5qZdYWSWt8amwn0R8TaiNgMLAXm1FwzB1ic7C8DZklSRLwSEXdQScqtNbvVC83MusIwJn1UL92QbD1V7zSFypyKbQaSYwx1TTJbeT0weSTN9gM8MyuWYYymqF66IWtOxmZWLO0bTbEO2LPq9dTk2FDXDEgaR2UxtRdHEqzzyTiLqYlbWi7TtFdGQ2pm9l2cSdz7jzk1k7gH3XBJ6jFj4y9TjwmgnXbNJC7jJ2QTtx3aN854FbC3pOlUku484ISaa/qA+cBdwFxgRUSMaHE094zNrFjK7ekARsRWSacBtwBl4PKIWCPpXGB1RPQBi4AlkvqBl6gkbAAkPQXsAoyXdAxwWEQ8Ui+ek7GZFUsb/0KNiOXA8ppjZ1XtbwKOrXPvtOHEcjI2s2LxdGgzsxzo0unQTsZmVizuGZuZ5UC5O9Nad7bazKyOVhYAyiMnYzMrFteMzcxywD1jM7MccM/YzCwHurRn3PBXiKQvSNqz0TVmZrlSLre+5Uiz/vw3gZWSfizpFElvS6NRZmYjNoz1jPOkWWvWUlk27pvA/sAjkm6WNF/SzvVuql6wufe6m9rYXDOzJqTWtxxpVjOOiBgEbgVulbQdle+DOh44Dxiyp1y9YPPg6ptGtJycmdnI5CvJtqpZMn7Tp4qILVTW7+yTNLFjrTIzG6mc9Xhb1SwZH1fvRES82ua2mJmNXhGTcUT8JK2GmJm1Rc4ezLXK44zNrFi6s2PsZGxmRdOd2djJ2MyKpYg1YzOzruNkbGaWA36AZ2aWB+4Zm5llz2UKM7MccDKuIzJYmqK8XfoxAQbfyCbuuPGZhN33ewsziZsFTXprNoG3bs4m7paM4raFk7GZWeb8haRmZnng0RRmZjngnrGZWQ44GZuZ5YGTsZlZ9twzNjPLge7MxU7GZlYwHk1hZpYDLlOYmeVBdybj7uzPm5nVI7W+NX0rzZb0uKR+SWcMcX6CpKuT8yslTas699Xk+OOSPtYsVsNkLGm8pBMlfTR5fYKkCyWdKimjBSDMzBpoUzKWVAYuAg4HZgDHS5pRc9kC4OWI2Au4AFiY3DsDmAfsA8wG/jl5v7qalSm+m1wzUdJ8YCfgOmAWMBOY3+R+M7N0te8B3kygPyLWAkhaCswBHqm6Zg5wdrK/DLhQlcUx5gBLI+J14ElJ/cn73VUvWLNk/N6I+ANJ44B1wB4R8YakfwMerHeTpB6gJ3l5UkT0NolT931Geu9ojKW4Y+mzjrW4Y+mzvsnESS0XjWtyFUBvVdunAE9XnRsADqx5i19fExFbJa0HJifH7665d0qjtjT7FVKSNB7YGZgITEqOTwDqlikiojciDki20fyj9DS/pCPGUtyx9FnHWtyx9FlHpCZXjTZfjUqznvEi4DGgDJwJXCNpLfB+YGmH22ZmlqV1wJ5Vr6cmx4a6ZiCpIEwCXmzx3jdp2DOOiAuAg4APRMQ/AZ8EbgEWRMQ5TT+KmVn3WgXsLWl6UiGYB/TVXNPHb56dzQVWREQkx+cloy2mA3sD9zQK1nSccUT8vGr/l1SK1GnJ6k+GsRR3LH3WsRZ3LH3WtktqwKdR6YCWgcsjYo2kc4HVEdFHpXqwJHlA9xKVhE1y3fepPOzbCpwaEQ2/CkiRxdcimZnZm3jSh5lZDjgZm5nlQG6TcbNpiB2Kebmk5yQ9nEa8JOaekm6X9IikNZJOTynu9pLukfRgEje1B7KSypLul3RjWjGTuE9JekjSA5JWpxRzV0nLJD0m6VFJH0gh5ruSz7ht2yDpi52Om8T+i+Tn6WFJV0naPo24RZDLmnEybfAnwKFUBkuvAo6PiEca3jj6uAcDG4ErIuI9nYxVFXN3YPeIuE/SzsC9wDEpfFYBO0bExmRq+x3A6RFxd5Nb2xH7S8ABwC4RcWSn41XFfQo4ICJeSDHmYuDHEXFZ8kR+YvIgPK34ZSpDqg6MiJ92ONYUKj9HMyLiteQB1vKI+NdOxi2KvPaMfz0NMSI2UxnTPKfTQSPiv6k8EU1NRDwTEfcl+78CHqXJTJ02xY2I2Ji83C7ZOv6bWdJU4Ajgsk7HypqkScDBVJ64ExGb00zEiVnAE51OxFXGATskY24nAj9vcr0l8pqMh5qG2PEElbVkxaf3AStTileW9ADwHHBbRKQR9x+BvwIGU4hVK4BbJd2bTIPttOnA88B3k7LMZZJ2TCFutXnAVWkEioh1wHnAz4BngPURcWsasYsgr8l4zJG0E3At8MWI2JBGzIh4IyL2pTI7aKakjpZmJB0JPBcR93YyTgMHRcR+VFbhOjUpS3XSOGA/4OKIeB/wCpDK8w+orLoIHA1ck1K8t1D5C3Y6sAewo6RPpRG7CPKajIc9lbCbJTXba4ErI+K6tOMnfzrfTmWpv076EHB0UrtdCnwkWXQqFUnPjYh4DrieSjmskwaAgaq/OJZRSc5pORy4LyKeTSneR4EnI+L5iNhCZYXHD6YUu+vlNRm3Mg2xEJIHaYuARyPi/BTjvk3Srsn+DlQelj7WyZgR8dWImBoR06j8m66IiFR6TpJ2TB6QkpQKDgM6OmomIn4BPC3pXcmhWbx5+cVOO56UShSJnwHvlzQx+bmeReUZiLUgl1+7VG8aYqfjSroKOAR4q6QB4BsRsajDYT8EfBp4KKnfAvx1RCzvcNzdgcXJ0/YS8P2ISHWoWcp2A66v5AjGAd+LiJtTiPt54MqkU7EW+GwKMbf9wjkUOCmNeAARsVLSMuA+KlOA76cgU6PTkMuhbWZmY01eyxRmZmOKk7GZWQ44GZuZ5YCTsZlZDjgZm5nlgJOxmVkOOBmbmeXA/wN9nnXtpcFfVQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ens_preds = ensemble(train_pool, val_pool, num_samples=10, iters=1000, lr=0.2)\n",
    "ens_preds = np.asarray(ens_preds)\n",
    "data = np.mean(ens_preds, axis=0)[:, 1] # average estimated data uncertainty\n",
    "knowledge = np.var(ens_preds, axis=0)[:, 0] # estimated knowledge uncertainty\n",
    "\n",
    "print(\"Knowledge uncertainty via ensemble:\")            \n",
    "sns.heatmap(knowledge.reshape([num_cat, num_cat]), cmap=\"Reds\")\n",
    "plt.show()\n",
    "\n",
    "print(\"Average predicted data uncertainty:\")            \n",
    "sns.heatmap(data.reshape([num_cat, num_cat]), cmap=\"Reds\", vmin=0, vmax=0.05)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The models perfectly detected knowledge uncertainty inside the heart (we can see no traces of the original heart border on the figure). This illustrates how by estimating knowledge uncertainty, we can detect anomalous inputs.\n",
    "\n",
    "In practice, training an ensemble of several CatBoost models can be too expensive. If you need a cheaper method - you can use a virtual ensemble, obtained from a single trained model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def virt_ensemble(train_pool, val_pool, num_samples=10, iters=1000, lr=0.2, seed=0):\n",
    "    ens_preds = []\n",
    "    model = CatBoostRegressor(iterations=iters, learning_rate=lr, loss_function='RMSEWithUncertainty', \n",
    "                              posterior_sampling=True, verbose=False, random_seed=seed)\n",
    "    model.fit(train_pool, eval_set=val_pool)\n",
    "    ens_preds = model.virtual_ensembles_predict(test, prediction_type='VirtEnsembles', virtual_ensembles_count=num_samples)\n",
    "    return ens_preds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, knowledge uncertainty can be obtained by measuring the variance of the mean across the obtained models. Similarly, we can also compute the average estimated data uncertainty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Knowledge uncertainty via virtual ensemble:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAcuUlEQVR4nO3df5TVd33n8edr7oQoJIKLPTYyrLAN6k60JxqKXbW2DWpg9TDpMVZw2xM0depZWBNz1gXM7npM0+7BdYPdmNhMAyZmTQZE3Uy7NrFtaNVq+GGSRiGSjiSVQU2UAFmDQgbe+8f94F4m997vHeb+GD7zeni+xzuf7/f7/nzvF/KaD5/v936vIgIzM8tHV6cPwMzMmsvBbmaWGQe7mVlmHOxmZplxsJuZZcbBbmaWGQe7mVlmuos2kPQqoA+Yk5oOAEMR8WgrD8zMzM5M3RG7pDXAICBgR1oE3C1pbesPz8zMxkv1Pnkq6THgooh4bkz7NGB3RCyosV8/0A9w602fvKT/fSubdsA2icTJDnWsDvVrLTdj1oT/cD+gFzX8cfo/i2ey/MtUNBVzEngZ8M9j2i9I66qKiAFgAICjR/zMAjOzNioK9muAv5X0T8D+1PYvgQuB1a08MDOzM+E7QgqCPSLulfQKYBGnXzzdGREnWn1wZmbj1a0sZ1fGpfCumIg4CTzQhmMxM5uwLud6cbCbmZ1NPBXjYDezzHR5KsbBbmZ58YjdwW5mmfEcu4PdzDJT8lSMg93M8uKpGAe7mWXGUzEOdjPLjEfsDnabkA4NjTyH2np1Hg442fl2Rwe7mWWm27nuYDezvHgqxsFuZpnp8vP6HexmlhffFeNgN7PMeCrGwW5mmfGI3cFuZpnxF2042M0sM56K8Tkws8x0qfGliKQlkvZKGpa0tsr6cyVtTuu3S5pXsW5dat8r6bKimpLmpxrDqea01L5B0sNpeUzS4cJzUPzWzMzOHl2o4aUeSSXgZmAp0AuskNQ7ZrOrgEMRcSGwAVif9u0FlgMXAUuAWySVCmquBzakWodSbSLiQxFxcURcDNwEfLH4HJiZZaSJI/ZFwHBE7IuI48Ag0Ddmmz7gjvR6K7BYklL7YEQci4jHgeFUr2rNtM+lqQap5uVVjmkFcHfhOSh8a2ZmZ5GSGl8KzAH2V/w8ktqqbhMRo8ARYHadfWu1zwYOpxpV+5L0cmA+cH/RgTvYzSwr45mKkdQvaVfF0t/p469jObA1Ik4UbXjGd8VIem9EfOZM9zcza4Xx3MceEQPAQI3VB4C5FT/3pLZq24xI6gZmAgcL9q3WfhCYJak7jdqr9bUcWNXA25rQiP1jtVZU/hYc2HT7BLowMxufrnEsBXYCC9LdKtMoB+vQmG2GgCvT6yuA+yMiUvvydNfMfGABsKNWzbTPtlSDVPOeU51IehXwYuCbjZyDuiN2SY/UWgW8tNZ+p/0WPHrk7H2ws5mddZr18aSIGJW0GrgPKAGbImK3pOuBXRExBGwE7pQ0DDxNOahJ220B9gCjwKpTUyjVaqYu1wCDkm4AHkq1T1lO+WJsQ3mqettJehK4jPKtN6etAr4RES8r7MHBnq9OfRmDP1nYep36s50xa8J/uJtf/NKGD/7dh57M8i9T0Rz7XwLnRcTDY1dI+ruWHJGZ2QT4jpCCYI+Iq+qse0/zD8fMbGKyHIKPk58VY2ZZkafqHOxmlhfHuoPdzDLjOXYHu5llxjMxDnYzy4y/zNrBbmaZcaw72M0sM/7OUwe7mWVGHrM72JvpxFe3Fm/UAqXfeGdH+p1KV6lO3HJdR/otfeD6jvRLV6kz/TbB1PlbWZuD3cyy4qkYB7uZZcZ3xTjYzSwzjnUHu5llZgpd+qnJwW5mWXGuO9jNLDO+3dHBbmaZKTnXHexmlhfnup9waWaZ0Tj+V1hLWiJpr6RhSWurrD9X0ua0frukeRXr1qX2vZIuK6opaX6qMZxqTqtY97uS9kjaLemuouN2sJtZVqTGl/p1VAJuBpYCvcAKSb1jNrsKOBQRFwIbgPVp315gOXARsAS4RVKpoOZ6YEOqdSjVRtICYB3wxoi4CLim6Bw42M0sK13jWAosAoYjYl9EHAcGgb4x2/QBd6TXW4HFKn83Xx8wGBHHIuJxYDjVq1oz7XNpqkGqeXl6/X7g5og4BBARTzVyDuqS9CpJiyWdN6Z9SdG+ZmbtpvEsUr+kXRVLf0WpOcD+ip9HUhvVtomIUeAIMLvOvrXaZwOHU42xfb0CeIWkf5D0QCPZW/fiqaQPAquAR4GNkq6OiHvS6j8B7i3qwMysnbrG8QmliBgABlp3NE3RDSwAfgvoAb4q6TURcbjWDkUj9vcDl0TE5anof5F0dVpX8+xV/hYc2HR744dvZjZB4xmxFzgAzK34uSe1Vd1GUjcwEzhYZ99a7QeBWanG2L5GgKGIeC5N6zxGOehrKgr2roj4KUBEPEE53JdKupE65yUiBiJiYUQs7H/fyoIuzMyaR1LDS4GdwIJ0t8o0yhdDh8ZsMwRcmV5fAdwfEZHal6e7ZuZTDuIdtWqmfbalGqSap2ZH/jfl7EXSSyhPzeyrd+BFwf6kpItP/ZBC/h3AS4DXFOxrZtZ2XWp8qSfNd68G7qM8Hb0lInZLul7SsrTZRmC2pGHgWmBt2nc3sAXYQ3nKelVEnKhVM9VaA1ybas1OtUnbHpS0h3L4fzgiDtY7dpV/UdRYKfUAoxHxoyrr3hgR/1CvOABHj9TuIDP+oo18+Ys22mT6zAn/pXp47ryGM+fi/U9k+Ze47sXTiBips6441M3M2qzLN3H7kQJmlpcG5s6z52A3s6w41x3sZpYZj9gd7GaWGee6g93MMjOeT57mysFuZlnpKrpBfQpwsJtZVuTbHR3sZpYXXzzNNNi/8a9e3ZF+R44d70i/8JGO9Pq7I3s70u/2X2n/0ywu6f/ttvcJMPKmN3Sk355vbO9Iv83gXM802M1s6vKI3cFuZplxrjvYzSwzJd8V42A3s7x4KsbBbmaZca472M0sMw52B7uZZUaeY3ewm1lefPG0+DtPzczOKlLjS3EtLZG0V9KwpLVV1p8raXNav13SvIp161L7XkmXFdVMX3C9PbVvTl92jaSVkn4s6eG0/EHRcTvYzSwrkhpeCuqUgJuBpUAvsEJS75jNrgIORcSFwAZgfdq3F1gOXAQsAW6RVCqouR7YkGodSrVP2RwRF6fltqJz4GA3s6w0ccS+CBiOiH0RcRwYBPrGbNMH3JFebwUWq/wbow8YjIhjEfE4MJzqVa2Z9rk01SDVvPxMz0FhsEtaJOnX0uteSddK+rdn2qGZWSs1a8QOzAH2V/w8ktqqbhMRo8ARYHadfWu1zwYOpxrV+nqnpEckbZU0t+jA6wa7pI8C/xP4tKT/BnwKmAGslXRdUXEzs3Ybz4hdUr+kXRVLf6ePv4q/AOZFxK8Cf83//xdCTUV3xVwBXAycC/wI6ImIZyR9AtgO/HG1ndLJ6Qe49aZP0v++lY2+ATOzCekqNX5XTEQMAAM1Vh8AKkfHPamt2jYjkrqBmcDBgn2rtR8EZknqTqP2X2wfEQcrtr8N+HjR+yoK9tGIOAEclfS9iHgmdfQzSSdr7XTayTp6JIoOwsysWZr4SIGdwAJJ8ymH7HLgPWO2GQKuBL5JeSB8f0SEpCHgLkk3Ai8DFgA7AFWrmfbZlmoMppr3pPdzQUT8MPW3DHi06MCLgv24pOkRcRS45FSjpJlAzWA3M+uYJt3HHhGjklYD9wElYFNE7JZ0PbArIoaAjcCdkoaBpykHNWm7LcAeYBRYlQbJVKuZulwDDEq6AXgo1Qb4oKRlqc7TwMqiYy8K9jdHxLF0oJVBfg7l3yhmZpNLE58pEBFfBr48pu2/Vrz+OfCuGvv+MVWmq6vVTO37KN81M7Z9HbBuPMddN9hPhXqV9p8APxlPR2Zm7eCnO/qRAmaWm5I/nuNgN7Os+CFgDnYzy42nYhzsZpYXj9gd7GaWG4/YHexmlhmP2B3sZpYX+a4YB7uZZcZTMXkG+xv2facj/Z646xMd6bdTTn53e0f6XfiFP21/p499u/19Aj3fuKkj/Z7N5AF7nsFuZlOYR+wOdjPLi293dLCbWW48Ynewm1lefFeMg93McuOpGAe7mWXGUzEOdjPLi5/H7mA3s9x4KsbBbmZ58cVT8Bkws7xIjS+FpbRE0l5Jw5LWVll/rqTNaf12SfMq1q1L7XslXVZUU9L8VGM41Zw2pq93SgpJC4uO28FuZllRlxpe6taRSsDNwFKgF1ghqXfMZlcBhyLiQmADsD7t2wssBy4ClgC3SCoV1FwPbEi1DqXap47lfOBqoKHneIw72CV9drz7mJm1TfNG7IuA4YjYFxHHgUGgb8w2fcAd6fVWYLHKV2/7gMGIOBYRjwPDqV7VmmmfS1MNUs3LK/r5I8rB//NGTkHdOXZJQ2ObgN+WNAsgIpY10omZWduM4+KppH6gv6JpICIG0us5wP6KdSPA68eU+MU2ETEq6QgwO7U/MGbfOel1tZqzgcMRMTp2e0mvA+ZGxP+R9OFG3lfRxdMeYA9wGxCUg30h8D/q7VR5sm696ZP0v29lI8diZjZh47ndMYX4QOGGHSKpC7gRWDme/YqCfSHleZ3rgA9HxMOSfhYRf19vp9NO1tEjMZ4DMjObkObdFXMAmFvxc09qq7bNiKRuYCZwsGDfau0HgVmSutOo/VT7+cCrgb9Lv7B+GRiStCwidtU68LpnICJORsQG4L3AdZI+hW+RNLPJrHlz7DuBBelulWmUL4aOnZ4eAq5Mr68A7o+ISO3L010z84EFwI5aNdM+21INUs17IuJIRLwkIuZFxDzK0zt1Qx0aDOmIGAHeJentwDON7GNm1hFN+uRpmjNfDdwHlIBNEbFb0vXArogYAjYCd0oaBp6mHNSk7bZQnsoeBVZFxIny4T2/ZupyDTAo6QbgoVT7jKj8i6KFptBUzFT7BiW99k0d6Td+/mz7O+3QNyiV3n1NR/rtmOkzJ5zKox/6nYYzp3vDl7L8mKqnVcwsL35WjIPdzDLjYHewm1lmSqVOH0HHOdjNLC8esTvYzSwzDnYHu5llxsHuYDezzHT5obUOdjPLi4O9DcHe6g9AVdOhf4qV3vMfO9LviR1/1ZF+4+APO9IvL3hh27vs2AeFOvHfD5zd0xln87E3iUfsZpYVecTuYDezzHjE7mA3s8w42B3sZpYZB7uD3cwy40cKONjNLDMesTvYzSwzDnYHu5llxrc7OtjNLDMesdf/Mmszs7NO877MGklLJO2VNCxpbZX150ranNZvlzSvYt261L5X0mVFNdMXXG9P7ZvTl10j6QOSvi3pYUlfl9RbdNwOdjPLS6nU+FKHpBJwM7AU6AVWVAnVq4BDEXEhsAFYn/btpfzF1hcBS4BbJJUKaq4HNqRah1JtgLsi4jURcTHwceDGolMwrmCX9CZJ10p623j2MzNrm+aN2BcBwxGxLyKOA4NA35ht+oA70uutwGJJSu2DEXEsIh4HhlO9qjXTPpemGqSalwNExDMV/c0ACh8gVDfYJe2oeP1+4FPA+cBHq/2zxMys45oX7HOA/RU/j6S2qttExChwBJhdZ99a7bOBw6nG8/qStErS9yiP2D9YdOBFI/ZzKl73A2+NiI8BbwP+Xa2dJPVL2iVp18Cm24uOwcysebq6Gl4qsyot/Z0+/Goi4uaI+BVgDfCfi7YvuiumS9KLKf8CUET8OHXyrKTRWjtFxAAwAMCzhzv03FEzm5LGcVfMaVn1fAeAuRU/96S2atuMSOoGZgIHC/at1n4QmCWpO43aq/UF5ambTxe8rcIR+0zgW8Au4F9IugBA0nmA7ykys8mnq9T4Ut9OYEG6W2Ua5YuhQ2O2GQKuTK+vAO6PiEjty9NdM/OBBcCOWjXTPttSDVLNewAkLajo7+3APxUdeN0Re0TMq7HqJPA7RcXNzNquqzljzogYlbQauA8oAZsiYrek64FdETEEbATulDQMPE05qEnbbQH2AKPAqog4AVCtZupyDTAo6QbgoVQbYLWktwDPUb5b5tQvkpoUrf6Glk5MxUyxDyh06huUOP7zzvTbiW9QWrik7X0CU+8blKbPnHDHJ279SMMnrfSHf5JlWPiTp2aWlyk2sKvGwW5mefGzYhzsZpYZj9gd7GaWmeK7XbLnYDezvHgqxsFuZpnxVIyD3cwyI4/YHexmlpcmfUDpbNb6YPc/i1qutGhpR/r9wIy5xRu1wJ89u794o1z4v5/x88VTj9jNLDOeinGwm1lmPBXjYDezzHj6ysFuZpnxVIyD3cwy46kYB7uZZcZ3xTjYzSwznopxsJtZZjwV42A3s8x4xO5gN7PM+HZH/KvNzPLS1dX4UkDSEkl7JQ1LWltl/bmSNqf12yXNq1i3LrXvlXRZUU1J81ON4VRzWmq/VtIeSY9I+ltJLy88BYXvzMzsbNJVanypQ1IJuBlYCvQCKyT1jtnsKuBQRFwIbADWp317geXARcAS4BZJpYKa64ENqdahVBvgIWBhRPwqsBX4eOEpKHhjr5f0ovT6hZI+JukvJK2XNLOouJlZ20mNL/UtAoYjYl9EHAcGgb4x2/QBd6TXW4HFkpTaByPiWEQ8DgynelVrpn0uTTVINS8HiIhtEXE0tT8A9BQdeNGIfRNwquCfAjMp/1Y5Cnym1k6S+iXtkrRrYNPtRcdgZtY845iKqcyqtPRXVJoDVD5KdCS1UW2biBgFjgCz6+xbq302cDjVqNUXlEfxf1V0CoounnZVdLQwIl6XXn9d0sO1doqIAWAAgKNHouggzMyaZhwXT0/LqklO0u8BC4HfLNq2aMT+HUnvTa//UdLC1MErgOcmdJRmZq2grsaX+g4AlV860JPaqm4jqZvyrMbBOvvWaj8IzEo1nteXpLcA1wHLIuJY0YEXvbM/AH5T0vcoT/R/U9I+4M/TOjOzyaVJF0+BncCCdLfKNMoXQ4fGbDMEXJleXwHcHxGR2penu2bmAwuAHbVqpn22pRqkmvcASHotcCvlUH+qkVNQdyomIo4AK9MF1Plp+5GIeLKR4mZmbdekT55GxKik1cB9QAnYFBG7JV0P7IqIIWAjcKekYeBpykFN2m4LsAcYBVZFxAmAajVTl2uAQUk3UL4TZmNq/+/AecDny9dY+X5ELKt37Cr/omghz7Fny1+NZ003feaEU/nE1z7fcOaUfuNdWX6ayZ88NbO8+JOnDnYzy4yfFeNgN7O8yCN2B7uZZabLseYzYGZ58fPYHexmlhnPsTvYzSwznmN3sJtZZjxid7DbmfMHhWxS8ojdwW5mmSkVPgMmew52M8uLp2Ic7GaWGU/FONjNLDMesTvYzSwzHrE72M0sMyXHms+AmWXFDwFzsJtZbjzH7mA3s8x4xF74ZdZmZmcXdTW+FJWSlkjaK2lY0toq68+VtDmt3y5pXsW6dal9r6TLimqmL7jento3py+7RtKbJT0oaVTSFTTAwW5meZEaX+qWUQm4GVgK9AIrJPWO2ewq4FBEXAhsANanfXspf7H1RcAS4BZJpYKa64ENqdahVBvg+8BK4K5GT0HdYJf0QUmd+cZiM7MzUSo1vtS3CBiOiH0RcRwYBPrGbNMH3JFebwUWq3z1tg8YjIhjEfE4MJzqVa2Z9rk01SDVvBwgIp6IiEeAk42egqIR+x8B2yV9TdK/l/RLjRY2M+uIcUzFSOqXtKti6a+oNAeofNLdSGqj2jYRMQocAWbX2bdW+2zgcKpRq6+GFQX7PqCHcsBfAuyRdK+kKyWdX2unypM1sOn2Mz02M7PxG8dUTEQMRMTCimWg04ffDEV3xUREnAS+AnxF0jmU54ZWAJ8Aqo7g08kpn6CjR6JpR2tmVqhpd8UcACqnontSW7VtRiR1AzOBgwX7Vms/CMyS1J1G7dX6aljRiP20MxQRz0XEUESsAF5+pp2ambVMky6eAjuBBelulWmUL4YOjdlmCLgyvb4CuD8iIrUvT3fNzAcWADtq1Uz7bEs1SDXvOdNTUBTs7661IiKOnmmnZmYt06RgTyPn1cB9wKPAlojYLel6ScvSZhuB2ZKGgWuBtWnf3cAWYA9wL7AqIk7UqplqrQGuTbVmp9pI+jVJI8C7gFslndq+9iko/6JoIU/FmFmjps+c8DxKHNjbcOZoziuz/DSTP3lqZnnJMqrHx8FuZplxsjvYzSwvflaMg93MMuNgd7CbWWb82F4Hu5nlxiN2B7uZ5cVTMQ52M8uMg93BbjaptfoDhLWc1eF4Nh97czjYzSwr/jJrB7uZ5cZ3xTjYzSwzHrE72M0sMw52B7uZ5cbB7mA3s7x4xO5gN7PMONcd7GaWGd8V42A3s8x4KsbBbma5cbA72M0sLx6x1w92SdOA5cAPIuJvJL0HeAPlb9ceiIjn2nCMZmaNc7CjqPOQIUmfoxz+04HDwHnAF4HFad8rC3s4eqRDTzEyy8BUewjY9JkT73g8mdOM/iajiKi5AI+k/+8GngRK6WedWldjv35gV1r66/VR0P8Z7zuRZSr1O5Xe61Trdyq9Vy+nL0Uj9u8ArwNmAN8HXh4RT0t6AfBQRPzr5vx6qdn/rohY2Mo+pnq/U+m9TrV+p9J7tdMVXTzdCHwXKAHXAZ+XtA/4dWCwxcdmZmZnoG6wR8QGSZvT6x9I+izwFuDPI2JHOw7QzMzGp/B2x4j4QcXrw8DWlh7R6Qba2NdU7Xcqvdep1u9Ueq9Woe4cu5mZnX38UAUzs8w42M3MMjNpg13SEkl7JQ1LWtumPjdJeird5tkWkuZK2iZpj6Tdkq5uU78vkLRD0j+mfj/Wjn5T3yVJD0n6y3b1mfp9QtK3JT0saVeb+pwlaauk70p6VNK/aUOfr0zv8dTyjKRrWt1v6vtD6e/TdyTdnW6NtjablHPskkrAY8BbgRFgJ7AiIva0uN83Az8FPhsRr25lXxV9XgBcEBEPSjof+BZweRveq4AZEfFTSecAXweujogHWtlv6vtaYCHwooh4R6v7q+j3CWBhRPykjX3eAXwtIm5Lj+iYnm5CaFf/JeAA8PqI+OcW9zWH8t+j3oj4maQtwJcj4vZW9mvPN1lH7IuA4YjYFxHHKd8z39fqTiPiq8DTre5nTJ8/jIgH0+v/S/k5PHPa0G9ExE/Tj+ekpeW/5SX1AG8Hbmt1X50maSbwZsqfByEijrcz1JPFwPdaHeoVuoEXSjr1KJIfFGxvLTBZg30OsL/i5xHaEHadJmke8Fpge5v6K0l6GHgK+OuIaEe/nwT+E3CyDX2NFcBXJH1LUn8b+psP/Bj4TJp6uk3SjDb0W2k5cHc7OoqIA8AnKH9K/YfAkYj4Sjv6ttNN1mCfciSdB3wBuCYinmlHnxFxIiIuBnqARZJaOv0k6R3AUxHxrVb2U8ebIuJ1wFJgVZp6a6Vuyo/k+HREvBZ4FmjL9SL4xdNZlwGfb1N/L6b8L+v5wMuAGZJ+rx192+kma7AfAOZW/NyT2rKU5ri/AHwuIr7Y7v7T9MA2YEmLu3ojsCzNdQ8Cl0r6Xy3u8xfSiJKIeAr4EuUpv1YaAUYq/iW0lXLQt8tS4MGIeLJN/b0FeDwifhzlR3p/kfJjvq3NJmuw7wQWSJpf8Uz4oQ4fU0uki5gbgUcj4sY29vtLkmal1y+kfKH6u63sMyLWRURPRMyj/Gd6f0S0ZUQnaUa6OE2aDnkb0NK7nyLiR8B+Sa9MTYuBll4UH2MFbZqGSb4P/Lqk6env9WLK14yszSblNyhFxKik1cB9lB9Atikidre6X0l3A78FvETSCPDRiNjY4m7fCPw+8O003w3wkYj4cov7vQC4I9010QVsiYi23n7YZi8FvlTOG7qBuyLi3jb0+x+Az6UByj7gvW3o89Qvr7cCf9iO/gAiYrukrcCDwCjwEH68QEdMytsdzczszE3WqRgzMztDDnYzs8w42M3MMuNgNzPLjIPdzCwzDnYzs8w42M3MMvP/AIueRvyU1Is+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average predicted data uncertainty:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAD8CAYAAACihcXDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAXDElEQVR4nO3df7RdZX3n8ffn3oQQQEMXWhdNkKQLxAa7Fj8ywRmV0UZomFpDK4zB6UhdjKnj0Oq4OlOmtajY6Rq6HJmuJTPtHcEiY/kVQG81FVR0AAdDEgiV8MMJEZsLivKjYORHCPczf5wd5/R67znn3pyz9z77fl5r7ZV99t7nfp8DWZ/75NnPfo5sExER1RqpugEREZEwjoiohYRxREQNJIwjImogYRwRUQMJ44iIGkgYR0TMQNJaSQ9K2inpgmnOL5J0TXF+s6TlxfHlkp6TtL3Y/qJbrQU9NOa1wDpgaXHoEWDc9v2z+VAREcNE0ihwKXAaMAFskTRu+762y84DnrJ9jKT1wMXAO4tzD9k+odd6HXvGkv4AuBoQcGexCbhqut8SERENshrYaXuX7b20snDdlGvWAVcU+xuBNZI0l2LdesbnAcfbfrH9oKRPAjuA/zLdmyRtADYA/MXH//DkDet/cy5tmzO/8Gyp9fbTQQdXUtfP7amkrhYfVkldv/Bc+UUXLCy/JqBFiyup62d/XEndkePfNKcga/c+vbznx4r/kh//DkVWFcZsjxX7S4HdbecmgFOm/IifXmN7n6SngSOKcysk3Q08A3zY9m2d2tItjCeBXwC+N+X4kcW5aRUfZgzAO7fleeuIqKX2rOqz7wOvtv2EpJOBz0s63vYzM72hWxh/EPiapP/L//8N8WrgGOD8frQ4IqKf+jgr4RHgqLbXy4pj010zIWkBsAR4wq1Ff14AsL1N0kPAa4CtMxXrGMa2vyzpNbTGTtpv4G2x/VLPHykioiQL5jZkO50twLGSVtDKvfXAu6ZcMw6cC9wBnAXcYtuSXgk8afslSb8IHAvs6tjubq2xPQl8a9YfIyKiAiN9yuJiDPh84CZgFLjc9g5JFwFbbY8DlwFXStoJPEkrsAFOBS6S9CKtId332X6yU72uYRwRMUz6+fCE7U3ApinHLmzbfx44e5r3XQ9cP5taCeOIaJSR/g1TlCphHBGNMqyPFSeMI6JR+jVmXLaEcUQ0ymiGKSIiqpdhioiIGsgwRUREDaRnPIPKFu1ZuKj0mv/9xDNKrwnw/m1frKRuJQv2AFQwJvjVN72j9JoAp91+QyV12be3mrp9kKltNVJFEEdEPSwYzixuZhhHxPyVYYqIiBoYYTi7xgnjiGiUzKaIiKiBDFNERNRAesYRETXQx8XlS5UwjohGyTBFREQNZJgiIqIGMrUtIqIG0jOOiKiB0YRxRET1hnWYYs43HiW9p58NiYjohxH1vtXJgcwC+dhMJyRtkLRV0tax68YPoERExOyMzGKrk47DFJL+bqZTwKtmep/tMWAMYHLHbZ5z6yIiZqlmHd6edRszfhXwq8BTU44L+D8DaVFExAFo6uLyXwQOs7196glJ3xhIiyIiDkDdhh961TGMbZ/X4dy7+t+ciIgDM5z94kxti4iGUUOHKSIihspwRnHCOCIappFjxhERw2ZIRykSxhHRLMP6OHTCOCIaZTijOGEcEQ1TtzUnepUwjohG0ZD2jQcexlqwcNAlfsbNbziz9JoAxy1eVEnd2089q5K6b7jlqkrqfunUf1l6zV//3o7SawLcfdyJldQ94dYbK6nbD8MZxekZR0TDDOswxbBOyYuImNYI6nnrRtJaSQ9K2inpgmnOL5J0TXF+s6TlU86/WtIeSb/fvd0REQ2iWWwdf440ClwKnAGsBM6RtHLKZecBT9k+BrgEuHjK+U8Cf9tLuxPGEdEoUu9bF6uBnbZ32d4LXA2sm3LNOuCKYn8jsEbF4hiSzgS+C/R0wyFhHBGNMpuecfu3EhXbhrYftRTY3fZ6ojjGdNfY3gc8DRwh6TDgD+jwjUhT5QZeRDTKbKa2tX8rUZ99FLjE9p5eV5FLGEdEo4z2bzbFI8BRba+XFcemu2ZC0gJgCfAEcApwlqQ/Aw4HJiU9b/tTMxVLGEdEo/RxZtsW4FhJK2iF7npg6pdqjAPnAncAZwG32Dbwpp+2R/oosKdTEEPCOCIapl9P4NneJ+l84CZgFLjc9g5JFwFbbY8DlwFXStoJPEkrsOckYRwRjdLPJTRtbwI2TTl2Ydv+88DZXX7GR3uplTCOiEYZ1iliXdst6bWS1hRTNdqPrx1csyIi5qZfD32UrWMYS/o94AvA7wL3Smqf8Pyng2xYRMRcjEg9b3XSrWf8XuBk22cCbwb+WNIHinMzfpL2idRj13y+Py2NiOjBsPaMu40Zj9jeA2D7YUlvBjZKOpoOn6V9IrUf/Jb71NaIiK56fciibrr1jB+TdML+F0Uwvw14BfDLg2xYRMRcjKj3rU669YzfDexrP1A8f/1uSX85sFZFRMyR6payPeoYxrYnOpz7Zv+bExFxYEaGdG5b5hlHRKMM65hxwjgiGmVIszhhHBHNkp5xREQNDGkWJ4wjolnq9mRdrxLGEdEoI02c2hYRMWyUqW0REdXLDbwZXPf6qd9sPXhn77i19JoAk3d9o5K6Vbnn9N+qpO7b7vhC+UX3vVh+TeDEB7ZVUvf6o19XSd13PPmDA/4ZQ5rF6RlHRLOkZxwRUQNDmsUJ44holtHMpoiIqF6GKSIiamBIszhhHBHNkjCOiKiBRi4uHxExbHIDLyKiBjJMERFRA5lNERFRA0Oaxd3DWNJqwLa3SFoJrAUesL1p4K2LiJilRvaMJX0EOANYIOkrwCnA14ELJJ1o+z+X0MaIiJ4NaRZ37RmfBZwALAJ+ACyz/YykTwCbgWnDWNIGYAPAexe/jLcuWty/FkdEdDAyOpxp3G0Z5n22X7L9LPCQ7WcAbD8HTM70JttjtlfZXpUgjogySep5q5NuPeO9kg4pwvjk/QclLaFDGEdEVKah84xPtf0CgO328F0InDuwVkVEzFXNery96hjG+4N4muOPA48PpEUREQegbsMPvco844holtHh/EbShHFENMqwLhQ0nL9CIiJmIvW+df1RWivpQUk7JV0wzflFkq4pzm+WtLw4vlrS9mK7R9JvdKuVnnFENEq/esaSRoFLgdOACWCLpHHb97Vddh7wlO1jJK0HLgbeCdwLrLK9T9KRwD2S/sb2vpnqpWccEc3Sv57xamCn7V229wJXA+umXLMOuKLY3wiskSTbz7YF78GAuxVLGEdEs4yo503SBklb27YNbT9pKbC77fVEcYzprinC92ngCABJp0jaAXwbeF+nXjFkmCIiGkazmE1hewwYG0Q7bG8Gjpf0S8AVkv7W9vMzXZ+ecUQ0S/+GKR4Bjmp7vaw4Nu01khYAS4An2i+wfT+wB3hdp2ID7xmfdfvGQZf4GXee9NbSawLc/A8/qaRuVT78nVsrqTt551dLrzn6L44rvSbA145eWUnd39wyvCvkqn9dzC3AsZJW0Ard9cC7plwzTutp5DtoLax2i20X79ld3MA7Gngt8HCnYhmmiIhm6dMTeEWQng/cBIwCl9veIekiYKvtceAy4EpJO4EnaQU2wBtpLTX8Iq11fN5fPLk8o4RxRDRKPx/6KL5EY9OUYxe27T8PnD3N+64ErpxNrYRxRDRL1qaIiKjebGZT1EnCOCKaZUjXpkgYR0SzZJgiIqJ6Wc84IqIOMkwREVG93MCLiKiDDFNERFRv3nzTh6TPDqIhERF90cdv+ihTx56xpPGph4C3SDocwPbbB9WwiIg5aWjPeBnwDPBJ4L8W24/b9qfVvmDz2HVT8zwiYnAk9bzVSbcx41XAB4A/Av6D7e2SnrP9vzu9qX3B5skdt3X9upGIiL5p4mwK25PAJZKuK/58rNt7IiIqVbMeb696ClbbE8DZkn6N1rBFREQ9NTmM97P9JeBLA2pLRMSBG2ngMEVExNCZDz3jiIjaSxhHRNTA6GjVLZiThHFENEt6xhERNZAwjoiogYRxREQNZGpbREQNJIxnoPL/w/yTr/116TUBnnvLOZXUvfpH1TwUOXnXNyqpy73bSy/5V+cdU3pNgHO3f6WSun7huUrq9mWAIcMUERHVU3rGERE1kJ5xREQNJIwjImogYRwRUQN5HDoiogbSM46IqIGEcUREDWRqW0REDaRnHBFRAwnjiIgaGNLZFLMaXJH0RkkfknT6oBoUEXFApN63rj9KayU9KGmnpAumOb9I0jXF+c2SlhfHT5O0TdK3iz9/pVutjmEs6c62/fcCnwJeBnxkuoZFRFSuT2EsaRS4FDgDWAmcI2nllMvOA56yfQxwCXBxcfxx4Ndt/zJwLnBlt2Z36xkvbNvfAJxm+2PA6cC/6vAhNkjaKmnr2LVf6NaGiIj+GRnpfetsNbDT9i7be4GrgXVTrlkHXFHsbwTWSJLtu20/WhzfASyWtKhTsW5jxiOSfo5WaMv2jwBs/0TSvpneZHsMGAOYvO+b7lIjIqJ/+ncDbymwu+31BHDKTNfY3ifpaeAIWj3j/d4B3GX7hU7FuoXxEmAbrWVGLelI29+XdBh9Wno0IqKvRnq/gSdpA61/9e83VnQm+0LS8bSGLrreZ+sYxraXz3BqEviNWbcsImLQRnrvJ7b/K34ajwBHtb1eVhyb7poJSQtodWCfAJC0DLgReLfth7o2u+dWt7H9rO3vzuW9EREDpZHet862AMdKWiHpIGA9MD7lmnFaN+gAzgJusW1JhwNfAi6w/c1emj2czw1GRMykT7MpbO8DzgduAu4HrrW9Q9JFkt5eXHYZcISkncCHgP2zzM4HjgEulLS92H6+U7089BERzdLHtSlsbwI2TTl2Ydv+88DZ07zvT4A/mU2thHFENEseh46IqIFZzKaok4RxRDRLltCMiKiBDFNERNRA9ylrtZQwjohmmcVDH3Uy+DCefGngJabSQQeXXhPgn99+QyV1l7zlZ2bWlGLk9WdUUnfLv/2z0mu+e/PUuf7l8AvPVVJXQ3oTDMgNvIiIWsgwRUREDWSYIiKiBjKbIiKiBjJMERFRAxmmiIiogcymiIiogQxTRETUQIYpIiJqID3jiIgayNS2iIgayBKaERE1MKSzKTr+CpF0iqSXF/uLJX1M0t9IuljSknKaGBExC336QtKydevPXw48W+z/ObAEuLg49pmZ3iRpg6StkraOXVfNalcRMU+NjPS+1Ui3YYqR4uuqAVbZPqnYv13S9pneZHsMGAOYvPdWH3gzIyJ6VLMeb6+6/Wq4V9J7iv17JK0CkPQa4MWBtiwiYi400vtWI916xv8G+HNJHwYeB+6QtBvYXZyLiKiXIb2B1zGMbT8N/HZxE29Fcf2E7cfKaFxExKw1+Qk8288A9wy4LRERB65mww+9yjzjiGiWIb2BlzCOiGZJzzgionpKzzgiogZGhjPWhrPVEREzafJsioiIoZEx44iIGsiYcUREDaRnPFOFhQMvMZX3VbRsRkWrQJ1wyzWV1N124lsqqbvqq58rv2hVf6f2Pl9N3SWvqKZuP6RnHBFRA6MNXJsiImLoDOkwxXC2OiJiJn38pg9JayU9KGmnpAumOb9I0jXF+c2SlhfHj5D0dUl7JH2ql2YnjCOiWfq0nrGkUeBS4AxgJXCOpJVTLjsPeMr2McAltL4JCeB54I+B3++12QnjiGiW/vWMVwM7be+yvRe4Glg35Zp1wBXF/kZgjSTZ/ont22mFck8SxhHRLKMLet7av6+z2Da0/aSltL5IY7+J4hjTXVN8Rd3TwBFzaXZu4EVEo8xmoaD27+usWsI4Ipqlf7MpHgGOanu9rDg23TUTkhYAS4An5lIswxQR0Sz9GzPeAhwraYWkg4D1wPiUa8aBc4v9s4BbbHsuzU7POCKapU89Y9v7JJ0P3ASMApfb3iHpImCr7XHgMuBKSTuBJ2kFdqsZ0sPAy4GDJJ0JnG77vpnqJYwjoln6+Di07U3ApinHLmzbfx44e4b3Lp9NrY6/QiT9nqSjOl0TEVEro6O9bzXSrT//cWCzpNskvV/SK8toVETEnPXpoY+ydWvNLlp3ED8OnAzcJ+nLks6V9LKZ3tQ+d2/s2s/3sbkREV308XHoMnUbM7btSeBm4GZJC2k9GngO8Alg2p5y+9y9yQfumNOdxYiIualXyPaqWxj/o09l+0VaUznGJR0ysFZFRMxVzXq8veoWxu+c6YTtZ/vcloiIA9fEMLb9nbIaEhHRFzW7MderzDOOiGYZzo5xwjgimmY40zhhHBHN0sQx44iIoZMwjoiogdzAi4iog/SMIyKql2GKiIgaSBjPYHJy4CV+xkhFY0ZVfFaABQsrKXvSjZdWUrcKWrS4krquKlj27a2mbl8kjCMiKjebLyStk4RxRDRLZlNERNRAesYRETWQMI6IqIOEcURE9dIzjoiogeHM4oRxRDRMZlNERNRAhikiIuogYRwRUb0m9owlHQSsBx61/VVJ7wL+GXA/MGb7xRLaGBHRuyaGMfCZ4ppDJJ0LHAbcAKwBVgPnDrZ5ERGzNKQ38LA94wb8XfHnAuAxYLR4rf3nZnjfBmBrsW3oVKNL/Tm/90C2+VR3Pn3W+VZ3Pn3WJmwq/uNNS9K9wEnAocDfA0fbflLSwcDdtn+pP78SZqy/1faqQdaY73Xn02edb3Xn02dtgm7DFJcBDwCjwB8B10naBbweuHrAbYuImDc6hrHtSyRdU+w/KumzwFuB/2n7zjIaGBExH3Sd2mb70bb9fwA2DrRF/9hYibXma9359FnnW9359FmHXscx44iIKMeQzgGJiGiWhHFERA3UNowlrZX0oKSdki4oqeblkn5YTOkrhaSjJH1d0n2Sdkj6QEl1D5Z0p6R7irofK6NuUXtU0t2SvlhWzaLuw5K+LWm7pK0l1Txc0kZJD0i6X9I/LaHmccVn3L89I+mDg65b1P73xd+neyVdVUyDjR7UcsxY0ijwHeA0YALYApxj+74B1z0V2AN81vbrBlmrreaRwJG275L0MmAbcGYJn1XAobb3SFoI3A58wPa3Blm3qP0hYBXwcttvG3S9troPA6tsP15izSuA22x/ulhe4JDiRnhZ9UeBR4BTbH9vwLWW0vp7tNL2c5KuBTbZ/qtB1m2KuvaMVwM7be+yvZfWnOZ1gy5q+1bgyUHXmVLz+7bvKvZ/TGvdj6Ul1LXtPcXLhcU28N/MkpYBvwZ8etC1qiZpCXAqrfn62N5bZhAX1gAPDTqI2ywAFktaABwCPNrl+ijUNYyXArvbXk9QQkBVTdJy4ERgc0n1RiVtB34IfMV2GXX/G/AfgckSak1l4GZJ2yRtKKHeCuBHwGeKYZlPSzq0hLrt1gNXlVHI9iPAJ2g9rft94GnbN5dRuwnqGsbzjqTDgOuBD9p+poyatl+yfQKwDFgtaaBDM5LeBvzQ9rZB1ungjbZPAs4A/l0xLDVIC2gtJ/A/bJ8I/AQo5f4H/HTVxbcD15VU7+do/Qt2BfALwKGSfquM2k1Q1zB+BDiq7fWy4lgjFWO21wOfs31D2fWLfzp/HVg74FJvAN5ejN1eDfyKpP814Jo/VfTcsP1D4EZaw2GDNAFMtP2LYyOtcC7LGcBdth8rqd5bge/a/pFby+veQGvJ3ehBXcN4C3CspBVtayqPV9ymgShupF0G3G/7kyXWfaWkw4v9xbRulj4wyJq2/5PtZbaX0/p/eovtUnpOkg4tbpBSDBWcDgx01oztHwC7JR1XHFoDDPTG7BTnUNIQReHvgddLOqT4e72G1j2Q6EEtv+nD9j5J5wM30Vqk6HLbOwZdV9JVwJuBV0iaAD5i+7IBl30D8K+BbxfjtwB/aHvTgOseCVxR3G0fAa61XepUs5K9CrixlREsAP7a9pdLqPu7wOeKTsUu4D0l1Nz/C+c04HfKqAdge7OkjcBdwD7gbvJodM9qObUtImK+qeswRUTEvJIwjoiogYRxREQNJIwjImogYRwRUQMJ44iIGkgYR0TUwP8D2PUXe3yi0wAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ens_preds = virt_ensemble(train_pool, val_pool, num_samples=10, iters=1000, lr=0.2)\n",
    "\n",
    "ens_preds = np.asarray(ens_preds)\n",
    "ens_preds = ens_preds\n",
    "\n",
    "data = np.mean(ens_preds, axis = 1)[:,1]\n",
    "knowledge = np.var(ens_preds, axis = 1)[:,0]\n",
    "\n",
    "print(\"Knowledge uncertainty via virtual ensemble:\")\n",
    "sns.heatmap(knowledge.reshape([num_cat,num_cat]), cmap=\"Reds\")\n",
    "plt.show()\n",
    "\n",
    "print(\"Average predicted data uncertainty:\")\n",
    "sns.heatmap(data.reshape([num_cat,num_cat]), cmap=\"Reds\", vmax = 0.05)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use prediction_type='TotalUncertainty' and get the same result easier. For this prediction type, CatBoost calculates and returns the following statistics using the virtual ensemble: [Mean Predictions, Knowledge Uncertainty, Data Uncertainty]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = CatBoostRegressor(iterations=1000, learning_rate=0.2, loss_function='RMSEWithUncertainty', \n",
    "                          posterior_sampling=True, verbose=False, random_seed=0)\n",
    "model.fit(train_pool, eval_set=val_pool)\n",
    "preds = model.virtual_ensembles_predict(test, prediction_type='TotalUncertainty', virtual_ensembles_count=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Knowledge uncertainty via virtual ensemble:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAcuUlEQVR4nO3df5TVd33n8edr7oQoJIKLPTYyrLAN6k60JxqKXbW2DWpg9TDpMVZw2xM0depZWBNz1gXM7npM0+7BdYPdmNhMAyZmTQZE3Uy7NrFtaNVq+GGSRiGSjiSVQU2UAFmDQgbe+8f94F4m997vHeb+GD7zeni+xzuf7/f7/nzvF/KaD5/v936vIgIzM8tHV6cPwMzMmsvBbmaWGQe7mVlmHOxmZplxsJuZZcbBbmaWGQe7mVlmuos2kPQqoA+Yk5oOAEMR8WgrD8zMzM5M3RG7pDXAICBgR1oE3C1pbesPz8zMxkv1Pnkq6THgooh4bkz7NGB3RCyosV8/0A9w602fvKT/fSubdsA2icTJDnWsDvVrLTdj1oT/cD+gFzX8cfo/i2ey/MtUNBVzEngZ8M9j2i9I66qKiAFgAICjR/zMAjOzNioK9muAv5X0T8D+1PYvgQuB1a08MDOzM+E7QgqCPSLulfQKYBGnXzzdGREnWn1wZmbj1a0sZ1fGpfCumIg4CTzQhmMxM5uwLud6cbCbmZ1NPBXjYDezzHR5KsbBbmZ58YjdwW5mmfEcu4PdzDJT8lSMg93M8uKpGAe7mWXGUzEOdjPLjEfsDnabkA4NjTyH2np1Hg442fl2Rwe7mWWm27nuYDezvHgqxsFuZpnp8vP6HexmlhffFeNgN7PMeCrGwW5mmfGI3cFuZpnxF2042M0sM56K8Tkws8x0qfGliKQlkvZKGpa0tsr6cyVtTuu3S5pXsW5dat8r6bKimpLmpxrDqea01L5B0sNpeUzS4cJzUPzWzMzOHl2o4aUeSSXgZmAp0AuskNQ7ZrOrgEMRcSGwAVif9u0FlgMXAUuAWySVCmquBzakWodSbSLiQxFxcURcDNwEfLH4HJiZZaSJI/ZFwHBE7IuI48Ag0Ddmmz7gjvR6K7BYklL7YEQci4jHgeFUr2rNtM+lqQap5uVVjmkFcHfhOSh8a2ZmZ5GSGl8KzAH2V/w8ktqqbhMRo8ARYHadfWu1zwYOpxpV+5L0cmA+cH/RgTvYzSwr45mKkdQvaVfF0t/p469jObA1Ik4UbXjGd8VIem9EfOZM9zcza4Xx3MceEQPAQI3VB4C5FT/3pLZq24xI6gZmAgcL9q3WfhCYJak7jdqr9bUcWNXA25rQiP1jtVZU/hYc2HT7BLowMxufrnEsBXYCC9LdKtMoB+vQmG2GgCvT6yuA+yMiUvvydNfMfGABsKNWzbTPtlSDVPOeU51IehXwYuCbjZyDuiN2SY/UWgW8tNZ+p/0WPHrk7H2ws5mddZr18aSIGJW0GrgPKAGbImK3pOuBXRExBGwE7pQ0DDxNOahJ220B9gCjwKpTUyjVaqYu1wCDkm4AHkq1T1lO+WJsQ3mqettJehK4jPKtN6etAr4RES8r7MHBnq9OfRmDP1nYep36s50xa8J/uJtf/NKGD/7dh57M8i9T0Rz7XwLnRcTDY1dI+ruWHJGZ2QT4jpCCYI+Iq+qse0/zD8fMbGKyHIKPk58VY2ZZkafqHOxmlhfHuoPdzDLjOXYHu5llxjMxDnYzy4y/zNrBbmaZcaw72M0sM/7OUwe7mWVGHrM72JvpxFe3Fm/UAqXfeGdH+p1KV6lO3HJdR/otfeD6jvRLV6kz/TbB1PlbWZuD3cyy4qkYB7uZZcZ3xTjYzSwzjnUHu5llZgpd+qnJwW5mWXGuO9jNLDO+3dHBbmaZKTnXHexmlhfnup9waWaZ0Tj+V1hLWiJpr6RhSWurrD9X0ua0frukeRXr1qX2vZIuK6opaX6qMZxqTqtY97uS9kjaLemuouN2sJtZVqTGl/p1VAJuBpYCvcAKSb1jNrsKOBQRFwIbgPVp315gOXARsAS4RVKpoOZ6YEOqdSjVRtICYB3wxoi4CLim6Bw42M0sK13jWAosAoYjYl9EHAcGgb4x2/QBd6TXW4HFKn83Xx8wGBHHIuJxYDjVq1oz7XNpqkGqeXl6/X7g5og4BBARTzVyDuqS9CpJiyWdN6Z9SdG+ZmbtpvEsUr+kXRVLf0WpOcD+ip9HUhvVtomIUeAIMLvOvrXaZwOHU42xfb0CeIWkf5D0QCPZW/fiqaQPAquAR4GNkq6OiHvS6j8B7i3qwMysnbrG8QmliBgABlp3NE3RDSwAfgvoAb4q6TURcbjWDkUj9vcDl0TE5anof5F0dVpX8+xV/hYc2HR744dvZjZB4xmxFzgAzK34uSe1Vd1GUjcwEzhYZ99a7QeBWanG2L5GgKGIeC5N6zxGOehrKgr2roj4KUBEPEE53JdKupE65yUiBiJiYUQs7H/fyoIuzMyaR1LDS4GdwIJ0t8o0yhdDh8ZsMwRcmV5fAdwfEZHal6e7ZuZTDuIdtWqmfbalGqSap2ZH/jfl7EXSSyhPzeyrd+BFwf6kpItP/ZBC/h3AS4DXFOxrZtZ2XWp8qSfNd68G7qM8Hb0lInZLul7SsrTZRmC2pGHgWmBt2nc3sAXYQ3nKelVEnKhVM9VaA1ybas1OtUnbHpS0h3L4fzgiDtY7dpV/UdRYKfUAoxHxoyrr3hgR/1CvOABHj9TuIDP+oo18+Ys22mT6zAn/pXp47ryGM+fi/U9k+Ze47sXTiBips6441M3M2qzLN3H7kQJmlpcG5s6z52A3s6w41x3sZpYZj9gd7GaWGee6g93MMjOeT57mysFuZlnpKrpBfQpwsJtZVuTbHR3sZpYXXzzNNNi/8a9e3ZF+R44d70i/8JGO9Pq7I3s70u/2X2n/0ywu6f/ttvcJMPKmN3Sk355vbO9Iv83gXM802M1s6vKI3cFuZplxrjvYzSwzJd8V42A3s7x4KsbBbmaZca472M0sMw52B7uZZUaeY3ewm1lefPG0+DtPzczOKlLjS3EtLZG0V9KwpLVV1p8raXNav13SvIp161L7XkmXFdVMX3C9PbVvTl92jaSVkn4s6eG0/EHRcTvYzSwrkhpeCuqUgJuBpUAvsEJS75jNrgIORcSFwAZgfdq3F1gOXAQsAW6RVCqouR7YkGodSrVP2RwRF6fltqJz4GA3s6w0ccS+CBiOiH0RcRwYBPrGbNMH3JFebwUWq/wbow8YjIhjEfE4MJzqVa2Z9rk01SDVvPxMz0FhsEtaJOnX0uteSddK+rdn2qGZWSs1a8QOzAH2V/w8ktqqbhMRo8ARYHadfWu1zwYOpxrV+nqnpEckbZU0t+jA6wa7pI8C/xP4tKT/BnwKmAGslXRdUXEzs3Ybz4hdUr+kXRVLf6ePv4q/AOZFxK8Cf83//xdCTUV3xVwBXAycC/wI6ImIZyR9AtgO/HG1ndLJ6Qe49aZP0v++lY2+ATOzCekqNX5XTEQMAAM1Vh8AKkfHPamt2jYjkrqBmcDBgn2rtR8EZknqTqP2X2wfEQcrtr8N+HjR+yoK9tGIOAEclfS9iHgmdfQzSSdr7XTayTp6JIoOwsysWZr4SIGdwAJJ8ymH7HLgPWO2GQKuBL5JeSB8f0SEpCHgLkk3Ai8DFgA7AFWrmfbZlmoMppr3pPdzQUT8MPW3DHi06MCLgv24pOkRcRS45FSjpJlAzWA3M+uYJt3HHhGjklYD9wElYFNE7JZ0PbArIoaAjcCdkoaBpykHNWm7LcAeYBRYlQbJVKuZulwDDEq6AXgo1Qb4oKRlqc7TwMqiYy8K9jdHxLF0oJVBfg7l3yhmZpNLE58pEBFfBr48pu2/Vrz+OfCuGvv+MVWmq6vVTO37KN81M7Z9HbBuPMddN9hPhXqV9p8APxlPR2Zm7eCnO/qRAmaWm5I/nuNgN7Os+CFgDnYzy42nYhzsZpYXj9gd7GaWG4/YHexmlhmP2B3sZpYX+a4YB7uZZcZTMXkG+xv2facj/Z646xMd6bdTTn53e0f6XfiFP21/p499u/19Aj3fuKkj/Z7N5AF7nsFuZlOYR+wOdjPLi293dLCbWW48Ynewm1lefFeMg93McuOpGAe7mWXGUzEOdjPLi5/H7mA3s9x4KsbBbmZ58cVT8Bkws7xIjS+FpbRE0l5Jw5LWVll/rqTNaf12SfMq1q1L7XslXVZUU9L8VGM41Zw2pq93SgpJC4uO28FuZllRlxpe6taRSsDNwFKgF1ghqXfMZlcBhyLiQmADsD7t2wssBy4ClgC3SCoV1FwPbEi1DqXap47lfOBqoKHneIw72CV9drz7mJm1TfNG7IuA4YjYFxHHgUGgb8w2fcAd6fVWYLHKV2/7gMGIOBYRjwPDqV7VmmmfS1MNUs3LK/r5I8rB//NGTkHdOXZJQ2ObgN+WNAsgIpY10omZWduM4+KppH6gv6JpICIG0us5wP6KdSPA68eU+MU2ETEq6QgwO7U/MGbfOel1tZqzgcMRMTp2e0mvA+ZGxP+R9OFG3lfRxdMeYA9wGxCUg30h8D/q7VR5sm696ZP0v29lI8diZjZh47ndMYX4QOGGHSKpC7gRWDme/YqCfSHleZ3rgA9HxMOSfhYRf19vp9NO1tEjMZ4DMjObkObdFXMAmFvxc09qq7bNiKRuYCZwsGDfau0HgVmSutOo/VT7+cCrgb9Lv7B+GRiStCwidtU68LpnICJORsQG4L3AdZI+hW+RNLPJrHlz7DuBBelulWmUL4aOnZ4eAq5Mr68A7o+ISO3L010z84EFwI5aNdM+21INUs17IuJIRLwkIuZFxDzK0zt1Qx0aDOmIGAHeJentwDON7GNm1hFN+uRpmjNfDdwHlIBNEbFb0vXArogYAjYCd0oaBp6mHNSk7bZQnsoeBVZFxIny4T2/ZupyDTAo6QbgoVT7jKj8i6KFptBUzFT7BiW99k0d6Td+/mz7O+3QNyiV3n1NR/rtmOkzJ5zKox/6nYYzp3vDl7L8mKqnVcwsL35WjIPdzDLjYHewm1lmSqVOH0HHOdjNLC8esTvYzSwzDnYHu5llxsHuYDezzHT5obUOdjPLi4O9DcHe6g9AVdOhf4qV3vMfO9LviR1/1ZF+4+APO9IvL3hh27vs2AeFOvHfD5zd0xln87E3iUfsZpYVecTuYDezzHjE7mA3s8w42B3sZpYZB7uD3cwy40cKONjNLDMesTvYzSwzDnYHu5llxrc7OtjNLDMesdf/Mmszs7NO877MGklLJO2VNCxpbZX150ranNZvlzSvYt261L5X0mVFNdMXXG9P7ZvTl10j6QOSvi3pYUlfl9RbdNwOdjPLS6nU+FKHpBJwM7AU6AVWVAnVq4BDEXEhsAFYn/btpfzF1hcBS4BbJJUKaq4HNqRah1JtgLsi4jURcTHwceDGolMwrmCX9CZJ10p623j2MzNrm+aN2BcBwxGxLyKOA4NA35ht+oA70uutwGJJSu2DEXEsIh4HhlO9qjXTPpemGqSalwNExDMV/c0ACh8gVDfYJe2oeP1+4FPA+cBHq/2zxMys45oX7HOA/RU/j6S2qttExChwBJhdZ99a7bOBw6nG8/qStErS9yiP2D9YdOBFI/ZzKl73A2+NiI8BbwP+Xa2dJPVL2iVp18Cm24uOwcysebq6Gl4qsyot/Z0+/Goi4uaI+BVgDfCfi7YvuiumS9KLKf8CUET8OHXyrKTRWjtFxAAwAMCzhzv03FEzm5LGcVfMaVn1fAeAuRU/96S2atuMSOoGZgIHC/at1n4QmCWpO43aq/UF5ambTxe8rcIR+0zgW8Au4F9IugBA0nmA7ykys8mnq9T4Ut9OYEG6W2Ua5YuhQ2O2GQKuTK+vAO6PiEjty9NdM/OBBcCOWjXTPttSDVLNewAkLajo7+3APxUdeN0Re0TMq7HqJPA7RcXNzNquqzljzogYlbQauA8oAZsiYrek64FdETEEbATulDQMPE05qEnbbQH2AKPAqog4AVCtZupyDTAo6QbgoVQbYLWktwDPUb5b5tQvkpoUrf6Glk5MxUyxDyh06huUOP7zzvTbiW9QWrik7X0CU+8blKbPnHDHJ279SMMnrfSHf5JlWPiTp2aWlyk2sKvGwW5mefGzYhzsZpYZj9gd7GaWmeK7XbLnYDezvHgqxsFuZpnxVIyD3cwyI4/YHexmlpcmfUDpbNb6YPc/i1qutGhpR/r9wIy5xRu1wJ89u794o1z4v5/x88VTj9jNLDOeinGwm1lmPBXjYDezzHj6ysFuZpnxVIyD3cwy46kYB7uZZcZ3xTjYzSwznopxsJtZZjwV42A3s8x4xO5gN7PM+HZH/KvNzPLS1dX4UkDSEkl7JQ1LWltl/bmSNqf12yXNq1i3LrXvlXRZUU1J81ON4VRzWmq/VtIeSY9I+ltJLy88BYXvzMzsbNJVanypQ1IJuBlYCvQCKyT1jtnsKuBQRFwIbADWp317geXARcAS4BZJpYKa64ENqdahVBvgIWBhRPwqsBX4eOEpKHhjr5f0ovT6hZI+JukvJK2XNLOouJlZ20mNL/UtAoYjYl9EHAcGgb4x2/QBd6TXW4HFkpTaByPiWEQ8DgynelVrpn0uTTVINS8HiIhtEXE0tT8A9BQdeNGIfRNwquCfAjMp/1Y5Cnym1k6S+iXtkrRrYNPtRcdgZtY845iKqcyqtPRXVJoDVD5KdCS1UW2biBgFjgCz6+xbq302cDjVqNUXlEfxf1V0CoounnZVdLQwIl6XXn9d0sO1doqIAWAAgKNHouggzMyaZhwXT0/LqklO0u8BC4HfLNq2aMT+HUnvTa//UdLC1MErgOcmdJRmZq2grsaX+g4AlV860JPaqm4jqZvyrMbBOvvWaj8IzEo1nteXpLcA1wHLIuJY0YEXvbM/AH5T0vcoT/R/U9I+4M/TOjOzyaVJF0+BncCCdLfKNMoXQ4fGbDMEXJleXwHcHxGR2penu2bmAwuAHbVqpn22pRqkmvcASHotcCvlUH+qkVNQdyomIo4AK9MF1Plp+5GIeLKR4mZmbdekT55GxKik1cB9QAnYFBG7JV0P7IqIIWAjcKekYeBpykFN2m4LsAcYBVZFxAmAajVTl2uAQUk3UL4TZmNq/+/AecDny9dY+X5ELKt37Cr/omghz7Fny1+NZ003feaEU/nE1z7fcOaUfuNdWX6ayZ88NbO8+JOnDnYzy4yfFeNgN7O8yCN2B7uZZabLseYzYGZ58fPYHexmlhnPsTvYzSwznmN3sJtZZjxid7DbmfMHhWxS8ojdwW5mmSkVPgMmew52M8uLp2Ic7GaWGU/FONjNLDMesTvYzSwzHrE72M0sMyXHms+AmWXFDwFzsJtZbjzH7mA3s8x4xF74ZdZmZmcXdTW+FJWSlkjaK2lY0toq68+VtDmt3y5pXsW6dal9r6TLimqmL7jento3py+7RtKbJT0oaVTSFTTAwW5meZEaX+qWUQm4GVgK9AIrJPWO2ewq4FBEXAhsANanfXspf7H1RcAS4BZJpYKa64ENqdahVBvg+8BK4K5GT0HdYJf0QUmd+cZiM7MzUSo1vtS3CBiOiH0RcRwYBPrGbNMH3JFebwUWq3z1tg8YjIhjEfE4MJzqVa2Z9rk01SDVvBwgIp6IiEeAk42egqIR+x8B2yV9TdK/l/RLjRY2M+uIcUzFSOqXtKti6a+oNAeofNLdSGqj2jYRMQocAWbX2bdW+2zgcKpRq6+GFQX7PqCHcsBfAuyRdK+kKyWdX2unypM1sOn2Mz02M7PxG8dUTEQMRMTCimWg04ffDEV3xUREnAS+AnxF0jmU54ZWAJ8Aqo7g08kpn6CjR6JpR2tmVqhpd8UcACqnontSW7VtRiR1AzOBgwX7Vms/CMyS1J1G7dX6aljRiP20MxQRz0XEUESsAF5+pp2ambVMky6eAjuBBelulWmUL4YOjdlmCLgyvb4CuD8iIrUvT3fNzAcWADtq1Uz7bEs1SDXvOdNTUBTs7661IiKOnmmnZmYt06RgTyPn1cB9wKPAlojYLel6ScvSZhuB2ZKGgWuBtWnf3cAWYA9wL7AqIk7UqplqrQGuTbVmp9pI+jVJI8C7gFslndq+9iko/6JoIU/FmFmjps+c8DxKHNjbcOZoziuz/DSTP3lqZnnJMqrHx8FuZplxsjvYzSwvflaMg93MMuNgd7CbWWb82F4Hu5nlxiN2B7uZ5cVTMQ52M8uMg93BbjaptfoDhLWc1eF4Nh97czjYzSwr/jJrB7uZ5cZ3xTjYzSwzHrE72M0sMw52B7uZ5cbB7mA3s7x4xO5gN7PMONcd7GaWGd8V42A3s8x4KsbBbma5cbA72M0sLx6x1w92SdOA5cAPIuJvJL0HeAPlb9ceiIjn2nCMZmaNc7CjqPOQIUmfoxz+04HDwHnAF4HFad8rC3s4eqRDTzEyy8BUewjY9JkT73g8mdOM/iajiKi5AI+k/+8GngRK6WedWldjv35gV1r66/VR0P8Z7zuRZSr1O5Xe61Trdyq9Vy+nL0Uj9u8ArwNmAN8HXh4RT0t6AfBQRPzr5vx6qdn/rohY2Mo+pnq/U+m9TrV+p9J7tdMVXTzdCHwXKAHXAZ+XtA/4dWCwxcdmZmZnoG6wR8QGSZvT6x9I+izwFuDPI2JHOw7QzMzGp/B2x4j4QcXrw8DWlh7R6Qba2NdU7Xcqvdep1u9Ueq9Woe4cu5mZnX38UAUzs8w42M3MMjNpg13SEkl7JQ1LWtumPjdJeird5tkWkuZK2iZpj6Tdkq5uU78vkLRD0j+mfj/Wjn5T3yVJD0n6y3b1mfp9QtK3JT0saVeb+pwlaauk70p6VNK/aUOfr0zv8dTyjKRrWt1v6vtD6e/TdyTdnW6NtjablHPskkrAY8BbgRFgJ7AiIva0uN83Az8FPhsRr25lXxV9XgBcEBEPSjof+BZweRveq4AZEfFTSecAXweujogHWtlv6vtaYCHwooh4R6v7q+j3CWBhRPykjX3eAXwtIm5Lj+iYnm5CaFf/JeAA8PqI+OcW9zWH8t+j3oj4maQtwJcj4vZW9mvPN1lH7IuA4YjYFxHHKd8z39fqTiPiq8DTre5nTJ8/jIgH0+v/S/k5PHPa0G9ExE/Tj+ekpeW/5SX1AG8Hbmt1X50maSbwZsqfByEijrcz1JPFwPdaHeoVuoEXSjr1KJIfFGxvLTBZg30OsL/i5xHaEHadJmke8Fpge5v6K0l6GHgK+OuIaEe/nwT+E3CyDX2NFcBXJH1LUn8b+psP/Bj4TJp6uk3SjDb0W2k5cHc7OoqIA8AnKH9K/YfAkYj4Sjv6ttNN1mCfciSdB3wBuCYinmlHnxFxIiIuBnqARZJaOv0k6R3AUxHxrVb2U8ebIuJ1wFJgVZp6a6Vuyo/k+HREvBZ4FmjL9SL4xdNZlwGfb1N/L6b8L+v5wMuAGZJ+rx192+kma7AfAOZW/NyT2rKU5ri/AHwuIr7Y7v7T9MA2YEmLu3ojsCzNdQ8Cl0r6Xy3u8xfSiJKIeAr4EuUpv1YaAUYq/iW0lXLQt8tS4MGIeLJN/b0FeDwifhzlR3p/kfJjvq3NJmuw7wQWSJpf8Uz4oQ4fU0uki5gbgUcj4sY29vtLkmal1y+kfKH6u63sMyLWRURPRMyj/Gd6f0S0ZUQnaUa6OE2aDnkb0NK7nyLiR8B+Sa9MTYuBll4UH2MFbZqGSb4P/Lqk6env9WLK14yszSblNyhFxKik1cB9lB9Atikidre6X0l3A78FvETSCPDRiNjY4m7fCPw+8O003w3wkYj4cov7vQC4I9010QVsiYi23n7YZi8FvlTOG7qBuyLi3jb0+x+Az6UByj7gvW3o89Qvr7cCf9iO/gAiYrukrcCDwCjwEH68QEdMytsdzczszE3WqRgzMztDDnYzs8w42M3MMuNgNzPLjIPdzCwzDnYzs8w42M3MMvP/AIueRvyU1Is+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average predicted data uncertainty:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAD8CAYAAACihcXDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAX/UlEQVR4nO3df7RdZX3n8ffn3pBIBEKHOi5MsMksQCfUDkpW0JE6TlM0VEpwCcvgmoqu6NUlWB3nR+nYItJOl8xiZJwFZcwysTG1Bow/GmgIOAUdaSEkQCiEH50QaUm0yK8JBgj5cT/zx9no8XrvOefenHP2Pvt+XmvtxT57P/t8n73I+p7nPvt5ni3bREREuYbKrkBERCQZR0RUQpJxREQFJBlHRFRAknFERAUkGUdEVECScUTEBCQtlfSIpB2SLhnn/CxJ1xXnN0uaXxyfL+lFSduK7X+1izWjg8q8HlgGzC0O7QY22H5oMjcVETFIJA0D1wBnAruALZI22H6wqdgK4FnbJ0paDlwBvLc496jtUzuN17JlLOn3gHWAgLuKTcDXxvuViIiokcXADts7be+nkQuXjSmzDFhT7K8HlkjSVIK1axmvAE6xfaD5oKTPA9uBz413kaQRYATg2k9/8rSR97xrKnWbugP7+hvvZTNmlhLWL+4tJa5mH1NKXF56sf8xh4f7HxNg1pHlxN33fClhhxafPaVE1uyjOqbjacVf5CcfochVhZW2Vxb7c4HHm87tAk4f8xU/LWP7oKQ9wHHFuQWS7gWeA/7A9vdb1aVdMh4FXgP8w5jjxxfnxlXczEqA0Xv/d+ZbR0QlNeeqLvsR8FrbT0s6Dfi2pFNsPzfRBe2S8SeBv5b0f/nZL8RrgROBi7tR44iIburiqITdwAlNn+cVx8Yrs0vSDGAO8LQbi/68BGD7bkmPAicDWycK1jIZ294k6WQafSfND/C22D7U8S1FRPTJjKl12Y5nC3CSpAU08t5y4H1jymwALgTuAM4DbrVtSa8CnrF9SNK/AE4Cdrasd7va2B4F7pz0bURElGCoS7m46AO+GLgZGAZW294u6XJgq+0NwCpgraQdwDM0EjbA24DLJR2g0aX7UdvPtIrXNhlHRAySbk6esL0R2Djm2KVN+/uA88e57hvANyYTK8k4ImplqHvdFH2VZBwRtTKo04qTjCOiVrrVZ9xvScYRUSvD6aaIiChfuikiIiog3RQRERWQlvFEylq0Z/iIvoe85tfHTs7pj4u+t7aUuKUs2FOSW876UClx37FpVSlxfWB/KXG7IUPbqqSERBwR1TBjMHNxTZNxRExb6aaIiKiAIQazaZxkHBG1ktEUEREVkG6KiIgKSMs4IqICuri4fF8lGUdEraSbIiKiAtJNERFRARnaFhFRAWkZR0RUwHCScURE+Qa1m2LKDx4lfbCbFYmI6IYhdb5VyeGMAvnsRCckjUjaKmnrym9tOowQERGTMzSJrUpadlNI+ruJTgGvnug62yuBlQCjd93oKdcuImKSKtbg7Vi7PuNXA+8Enh1zXMDf9qRGERGHoa6Ly98IHGV729gTkr7bkxpFRByGqnU/dKplMra9osW5ct4xFBHRwmC2izO0LSJqRjXtpoiIGCiDmYqTjCOiZmrZZxwRMWgGtJciyTgi6mVQp0MnGUdErQxmKk4yjoiaqdqaE51KMo6IWtGAto17n4yHj+h5iLE2vfMDfY8JcPLsmaXE/d5vfaiUuP/mxi+WEveGs/p/v+c8fFffYwJs/bUzSol72qYvlxK3G7qZiiUtBb4ADANfsv25MednAV8BTgOeBt5r+7Gm868FHgQus31lq1iDOgokImJc3VpCU9IwcA1wFrAQuEDSwjHFVgDP2j4RuAq4Ysz5zwM3dVTvTgpFRAyKIdTx1sZiYIftnbb3A+uAZWPKLAPWFPvrgSUqpgBKOhf4AbC9s3pHRNSIJrM1rb1ebCNNXzUXeLzp867iGOOVsX0Q2AMcJ+ko4Pdose77WHmAFxG1MplJH81rr3fZZcBVtvd2ulZGknFE1EoXH+DtBk5o+jyvODZemV2SZgBzaDzIOx04T9J/A44FRiXts331RMGSjCOiVro4tG0LcJKkBTSS7nJg7NLBG4ALgTuA84BbbRv49Z/WR7oM2NsqEUOScUTUzHCXcrHtg5IuBm6mMbRtte3tki4HttreAKwC1kraATxDI2FPSZJxRNRKN8cZ294IbBxz7NKm/X3A+W2+47JOYiUZR0StZAZeREQFZAnNiIgKGNTJE23rLen1kpYUg5ibjy/tXbUiIqZmMpM+qqRlMpb0u8BfAh8HHpDUPBXwT3pZsYiIqRiSOt6qpF3L+MPAabbPBd4O/KGkTxTnJryT5imGK7+5caJiERFdN6gt43Z9xkO29wLYfkzS24H1kn6FFvfSPMVw9O6b3aW6RkS01en046pp1zJ+QtKpL38oEvPZwC8Db+hlxSIipqJbS2j2W7uW8fuBg80HipWJ3i+pnJXFIyJaUNWybIdaJmPbu1qc+5vuVyci4vAMDejYtowzjohaGdQ+4yTjiKiVAc3FScYRUS9pGUdEVMCA5uIk44iol6rNrOtUknFE1MpQHYe2RUQMGmVoW0RE+fIAbwLXLfmdXof4Bcu33tT3mACjW28rJW5Zti77WClxf/s7X+l/0EP7+x8TWPR33y8l7vqTF5cS9/xnnzjs7xjQXJyWcUTUS1rGEREVMKC5OMk4IuplOKMpIiLKl26KiIgKGNBcnGQcEfWSZBwRUQG1XFw+ImLQ5AFeREQFpJsiIqICMpoiIqICBjQXt0/GkhYDtr1F0kJgKfCw7Y09r11ExCTVsmUs6TPAWcAMSd8BTgduAy6R9Ebb/7UPdYyI6NiA5uK2LePzgFOBWcA/AfNsPyfpSmAzMG4yljQCjACsOPIolsw8sns1johoYWh4MLNxu2WYD9o+ZPsF4FHbzwHYfhEYnegi2yttL7K9KIk4IvpJUsdblbRLxvslzS72T3v5oKQ5tEjGERGlGVLnWxuSlkp6RNIOSZeMc36WpOuK85slzS+OL5a0rdjuk/TuttVuc/5tRasY283J9wjgwrZ3EhHRb1LnW8uv0TBwDY3nZguBC4pBDM1WAM/aPhG4CriiOP4AsMj2qTQGPXxRUstu4ZbJ2PZLExx/yvb9Le8kIqIEXeymWAzssL3T9n5gHbBsTJllwJpifz2wRJJsv2D7YHH8FYDbBRvQV/dFRExgeKjjTdKIpK1N20jTN80FHm/6vKs4xnhliuS7BzgOQNLpkrYD9wMfbUrO48qkj4iolcksFGR7JbCyF/WwvRk4RdK/BNZIusn2vonKp2UcEfXSpT5jYDdwQtPnecWxccsUfcJzgKebC9h+CNgL/GqrYEnGEVErGlLHWxtbgJMkLZA0E1gObBhTZgM/G8xwHnCrbRfXzACQ9CvA64HHWgVLN0VE1EuXxg/bPijpYuBmYBhYbXu7pMuBrbY3AKuAtZJ2AM/QSNgAZ9CYqXyAxjDgj9l+qlW8JOOIqJcurmdcrMGzccyxS5v29wHnj3PdWmDtZGIlGUdErWh4MHtfk4wjol4qNs25Uz1Pxu/dtKrXIX7B37713L7HBLh1z/OlxC3Lp+/+y1Lieut3+x5zaH7LB+E9c8vrFpUS9z3f/YtS4naDBrNhnJZxRNRMWsYREeXL26EjIqogLeOIiPJlNEVERBWkmyIiogLSTRERUb6qvU6pU0nGEVEv6aaIiChfHuBFRFRBuikiIso3qJM+Jt2el/SVXlQkIqIruvemj75q2TKWNHZVewH/VtKxALbP6VXFIiKmpKYt43nAc8Dngf9ebD9p2h9X8xtXV357U7fqGhHRlqSOtypp12e8CPgE8GngP9neJulF299rdVHzG1dH77zBXalpREQn6jiawvYocJWkrxf/faLdNRERpapYi7dTHSVW27uA8yW9i0a3RURENdU5Gb/M9l8Bf9WjukREHL6hGnZTREQMnOnQMo6IqLwk44iIChgeLrsGU5JkHBH1kpZxREQFJBlHRFRAknFERAVkaFtERAUkGU+ghBWU3nLDtX2PCfDS2R8pJe71T/6klLjedns5cbff3/eYqz/+r/oeE+ADt19fSlxeerGcuN2QboqIiPIpLeOIiApIyzgiogKSjCMiKiDJOCKiAgZ0OvRg9nRHREykiy8klbRU0iOSdki6ZJzzsyRdV5zfLGl+cfxMSXdLur/472+0i5VkHBH10qVkLGkYuAY4C1gIXCBp4ZhiK4BnbZ8IXAVcURx/Cvht228ALgTWtqt2knFE1MvQUOdba4uBHbZ32t4PrAOWjSmzDFhT7K8HlkiS7Xtt/7A4vh04UtKsltWe1E1GRFTdJFrGzW+yL7aRpm+aCzze9HlXcYzxytg+COwBjhtT5j3APbZfalXtPMCLiHqZxGiK5jfZ96YqOoVG18U72pVNMo6IeuneaIrdwAlNn+cVx8Yrs0vSDGAO8DSApHnAt4D32360XbBJdVNIOkPSpyS1zfIREaXo3miKLcBJkhZImgksBzaMKbOBxgM6gPOAW21b0rE0Xt58ie2/6aTaLZOxpLua9j8MXA0cDXxmvGEeERGl61IyLvqALwZuBh4Crre9XdLlks4piq0CjpO0A/gU8HJevBg4EbhU0rZi++et4rXrpjiiaX8EONP2k5KuBO4EPjfeRUUn+AjAtZdcxMi7l7YJExHRJV1cKMj2RmDjmGOXNu3vA84f57o/Bv54MrHaJeMhSb9EowUt208WgZ6XdHCii5o7xUfvutGTqVBExGGp6XToOcDdgABLOt72jyQdVRyLiKiWocGcDt0yGdueP8GpUeDdXa9NRMThKuGFFt0wpaFttl8AftDlukREHD4N5ly2jDOOiHqpaZ9xRMRgyWuXIiIqIC3jiIgKqONoioiIgZNuioiICkg3RUREBWRoW0REBUynSR+TMlrC0hQzW77dpGfeftPqUuIeffaHS4k79NbfKiXuHf/hT/se8wN/vaZ9oV546cVy4g7oQzBgYOuelnFE1Eu6KSIiKiDdFBERFZDRFBERFZBuioiICkg3RUREBWQ0RUREBaSbIiKiAtJNERFRAWkZR0RUQIa2RURUQJbQjIiogAEdTdHyJ0TS6ZKOKfaPlPRZSTdIukLSnP5UMSJiEqTOtwpp155fDbxQ7H8BmANcURz78kQXSRqRtFXS1pXf3tSVikZEdGRoqPOtQtp1UwzZPljsL7L9pmL/dknbJrrI9kpgJcDonTeUsIZmRExbFWvxdqrdT8MDkj5Y7N8naRGApJOBAz2tWUTEVGio861C2rWMPwR8QdIfAE8Bd0h6HHi8OBcRUS0D+gCvZTK2vQf4QPEQb0FRfpftJ/pRuYiISavzDDzbzwH39bguERGHr2LdD53KOOOIqJcBfYCXZBwR9TKgLePBrHVExAQkdbx18F1LJT0iaYekS8Y5P0vSdcX5zZLmF8ePk3SbpL2Sru6k3knGEVEvQzM631qQNAxcA5wFLAQukLRwTLEVwLO2TwSuojEpDmAf8IfAf+y42p0WjIgYCEPqfGttMbDD9k7b+4F1wLIxZZYBa4r99cASSbL9vO3baSTlzqrdacGIiIEwiUkfzUs3FNtI0zfNpTGn4mW7imOMV6aYrbwHOG4q1c4DvIiol0mMpmheuqFsScYRUS/dG02xGzih6fO84th4ZXZJmkFjMbWnpxKs98l4Rgn5/uD+/seE0obUnHbDF0uJe9ebzy4l7ps3XNv/oIcO9T8m4AMddzl2lY75Z6XE7YrujTPeApwkaQGNpLsceN+YMhuAC4E7gPOAW21PaXG0tIwjol6Gu7M2he2Dki4GbgaGgdW2t0u6HNhqewOwClgraQfwDI2EDYCkx4BjgJmSzgXeYfvBieIlGUdEvXTxL1TbG4GNY45d2rS/Dzh/gmvnTyZWknFE1EumQ0dEVMCATodOMo6IeknLOCKiAoYHM60NZq0jIibQyQJAVZRkHBH1kj7jiIgKSMs4IqIC0jKOiKiAAW0Zt/wJkfS7kk5oVSYiolKGhzvfKqRde/6PgM2Svi/pY5Je1Y9KRURM2STWM66SdrXZSWPZuD8CTgMelLRJ0oWSjp7oouYFm1d+86YuVjciog2p861C2vUZ2/YocAtwi6QjaLwP6gLgSmDclnLzgs2jW2+a0nJyERFTU60k26l2yfjn7sr2ARrrd26QNLtntYqImKqKtXg71S4Zv3eiE7Zf6HJdIiIOXx2Tse2/71dFIiK6omIP5jqVccYRUS+D2TBOMo6IuhnMbJxkHBH1Usc+44iIgZNkHBFRAXmAFxFRBWkZR0SUL90UEREVkGQ8gdFDPQ/xC4ZKWhqvjHsFGDqilLCLvvonpcRltITlTmbO6n9MSvyD+8D+siJ3QZJxRETp8kLSiIgqyGiKiIgKSMs4IqICkowjIqogyTgionxpGUdEVMBg5uIk44iomYymiIiogHRTRERUwWAm48Fsz0dETETqfGv7VVoq6RFJOyRdMs75WZKuK85vljS/6dzvF8cfkfTOdrFaJmNJMyW9X9JvFp/fJ+lqSRdJKmdBhIiIVrqUjCUNA9cAZwELgQskLRxTbAXwrO0TgauAK4prFwLLgVOApcCfFt83oXbdFF8uysyWdCFwFPBNYAmwGLiwzfUREf3VvQd4i4EdtncCSFoHLAMebCqzDLis2F8PXK3G4hjLgHW2XwJ+IGlH8X13TBSsXTJ+g+1fkzQD2A28xvYhSX8O3DfRRZJGgJHi40dsr2wTZ8Lvmeq1h2M6xZ1O9zrd4k6ne/05s+d03Gk8JlcBrGyq+1zg8aZzu4DTx3zFT8vYPihpD3BccfzOMdfObVWXdj8hQ5JmAkcDs4E5xfFZwITdFLZX2l5UbIfzP2WkfZGemE5xp9O9Tre40+lep2RMrjrcfHVY2rWMVwEPA8PAp4GvS9oJvBlY1+O6RUSUaTdwQtPnecWx8crsKnoQ5gBPd3jtz2nZMrZ9FXAG8Bbb/xN4D3AzsML2Z9veSkTE4NoCnCRpQdFDsBzYMKbMBn727Ow84FbbLo4vL0ZbLABOAu5qFaztOGPbP2za/380Oqn7paw/GaZT3Ol0r9Mt7nS6164r+oAvptEAHQZW294u6XJgq+0NNHoP1hYP6J6hkbApyl1P42HfQeAi2y1fBaRGEo+IiDJl0kdERAUkGUdEVEBlk3G7aYg9irla0o8lPdCPeEXMEyTdJulBSdslfaJPcV8h6S5J9xVx+/ZAVtKwpHsl3divmEXcxyTdL2mbpK19inmspPWSHpb0kKS39CHm64p7fHl7TtInex23iP3vi39PD0j6mqRX9CNuHVSyz7iYNvj3wJk0BktvAS6w/WDLCw8/7tuAvcBXbP9qL2M1xTweON72PZKOBu4Gzu3DvQp4pe29xdT224FP2L6zzaXdiP0pYBFwjO2zex2vKe5jwCLbT/Ux5hrg+7a/VDyRn108CO9X/GEaQ6pOt/0PPY41l8a/o4W2XyweYG20/We9jFsXVW0Z/3Qaou39NMY0L+t1UNv/h8YT0b6x/SPb9xT7PwEeos1MnS7Fte29xccjiq3nv8yS5gHvAr7U61hlkzQHeBuNJ+7Y3t/PRFxYAjza60TcZAZwZDHmdjbwwzblo1DVZDzeNMSeJ6iyFSs+vRHY3Kd4w5K2AT8GvmO7H3H/B/CfgdE+xBrLwC2S7i6mwfbaAuBJ4MtFt8yXJL2yD3GbLQe+1o9AtncDVwL/CPwI2GP7ln7EroOqJuNpR9JRwDeAT9p+rh8xbR+yfSqN2UGLJfW0a0bS2cCPbd/dyzgtnGH7TTRW4bqo6JbqpRnAm4Brbb8ReB7oy/MPaKy6CJwDfL1P8X6Jxl+wC4DXAK+U9O/6EbsOqpqMJz2VcJAVfbbfAL5q+5v9jl/86XwbjaX+eumtwDlF3+064DeKRaf6omi5YfvHwLdodIf10i5gV9NfHOtpJOd+OQu4x/YTfYr3m8APbD9p+wCNFR7/dZ9iD7yqJuNOpiHWQvEgbRXwkO3P9zHuqyQdW+wfSeNh6cO9jGn7923Psz2fxv/TW233peUk6ZXFA1KKroJ3AD0dNWP7n4DHJb2uOLSEn19+sdcuoE9dFIV/BN4saXbx73oJjWcg0YFKvnZpommIvY4r6WvA24FflrQL+IztVT0O+1bgd4D7i/5bgP9ie2OP4x4PrCmetg8B19vu61CzPns18K1GjmAG8Be2N/Uh7seBrxaNip3AB/sQ8+UfnDOBj/QjHoDtzZLWA/fQmAJ8LzWZGt0PlRzaFhEx3VS1myIiYlpJMo6IqIAk44iICkgyjoiogCTjiIgKSDKOiKiAJOOIiAr4/wyIiWkQP0mKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mean_preds = preds[:,0] # mean values predicted by a virtual ensemble\n",
    "\n",
    "knowledge = preds[:,1] # knowledge uncertainty\n",
    "print(\"Knowledge uncertainty via virtual ensemble:\")\n",
    "sns.heatmap(knowledge.reshape([num_cat,num_cat]), cmap=\"Reds\")\n",
    "plt.show()\n",
    "\n",
    "data = preds[:,2] # estimated data uncertainty\n",
    "print(\"Average predicted data uncertainty:\")\n",
    "sns.heatmap(data.reshape([num_cat,num_cat]), cmap=\"Reds\", vmin = 0, vmax = 0.05)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
